text:
# Functions for visualizing graphs
export draw, drawPdf, drawPng

####################################################
# draw methods
####################################################

function graphviz(dot_graph::AbstractString; external_viewer::Symbol=:None)
    # Show a DOT graph
    validateGraphVizInstalled() # Show an error if GraphViz is not installed correctly
    if external_viewer == :iterm
        viewDotIniTerm(dot_graph)
    elseif external_viewer == :default
        viewDotExternal(dot_graph)
    else
        try
            # For iJulia notebook
            #display("text/html", "")
            display("text/html", dot2svg(dot_graph))
        catch
            viewDotExternal(dot_graph)
        end
    end
end

draw(factor_graph::FactorGraph; schedule=ScheduleEntry[], args...) = graphviz(genDot(nodes(factor_graph), edges(factor_graph), schedule=schedule); args...)
draw(; args...) = draw(currentGraph(); args...)

function drawPng(factor_graph::FactorGraph, filename::AbstractString)
    dot_graph = genDot(nodes(factor_graph), edges(factor_graph))
    png_graph = dot2png(dot_graph)
    png_file = open(filename, "w")
    write(png_file, png_graph)
end
drawPng(filename::AbstractString) = drawPng(currentGraph(), filename)

draw(nodes::Set{FactorNode}; args...) = graphviz(genDot(nodes, edges(nodes)); args...)
draw(nodes::Vector{FactorNode}; args...) = draw(Set(nodes); args...)

####################################################
# writePdf methods
####################################################

function dot2pdf(dot_graph::AbstractString, filename::AbstractString)
    # Generates a DOT graph and writes it to a pdf file
    validateGraphVizInstalled() # Show an error if GraphViz is not installed correctly
    open(`dot -Tpdf -o$(filename)`, "w", STDOUT) do io
        println(io, dot_graph)
    end
end

drawPdf(factor_graph::FactorGraph, filename::AbstractString) = dot2pdf(genDot(nodes(factor_graph), edges(factor_graph)), filename)
drawPdf(filename::AbstractString) = drawPdf(currentGraph(), filename)

drawPdf(nodes::Set{FactorNode}, filename::AbstractString) = dot2pdf(genDot(nodes, edges(nodes)), filename)
drawPdf(nodes::Vector{FactorNode}, filename::AbstractString) = drawPdf(Set(nodes), filename)

####################################################
# Internal functions
####################################################

function genDot(nodeset::Set{FactorNode}, edgeset::Set{Edge}; schedule::Schedule=ScheduleEntry[], external_edges::Set{Edge}=Set{Edge}())
    # Return a string representing the graph in DOT format
    # http://en.wikipedia.org/wiki/DOT_(graph_description_language)

    dot = "graph G{splines=true;sep=\"+25,25\";overlap=scalexy;nodesep=1.6;compound=true;\n"
    dot *= "\tnode [shape=box, width=1.0, height=1.0, fontsize=9];\n"
    dot *= "\tedge [fontsize=8, arrowhead=onormal];\n"

    # Draw nodes
    for node in nodeset
        if isDeltaConstraint(node)
            dot *= "\t$(objectid(node)) [label=\"$(node.id)\", style=filled, width=0.75, height=0.75]\n"
        elseif isa(node, Terminal)
            dot *= "\t$(objectid(node)) [label=\"Terminal $(node.id)\", style=filled, width=0.75, height=0.75]\n"
        else
            dot *= "\t$(objectid(node)) [label=\"$(slug(typeof(node)))\\n$(node.id)\"]\n"
        end
    end

    # Build dictionary for message labels
    msg_labels = Dict{Interface, String}()
    for (i, entry) in enumerate(condense(schedule))
        if entry.message_update_rule <: SumProductRule
            str = "($i)"
        elseif entry.message_update_rule <: Union{NaiveVariationalRule, StructuredVariationalRule}
            str = "(($i))"
        elseif entry.message_update_rule <: ExpectationPropagationRule
            str = "[$i]"
        else
            str = "?$i?"
        end
        msg_labels[entry.interface] = str
    end

    for edge in edgeset
        dot *= edgeDot(edge, msg_labels=msg_labels)
    end

    # Draw external edges
    if length(external_edges) > 0
        # Add invisible nodes to external edges
        external_nodes = setdiff(nodes(external_edges), nodeset)
        for node in external_nodes
            dot *= "\t$(objectid(node)) [label=\"\", shape=none,  width=0.15, height=0.15]\n"
        end

        # Add the external edges themselves
        for edge in external_edges
            dot *= edgeDot(edge, is_external_edge=true)
        end
    end

    dot *= "}";

    return dot
end

function edgeDot(edge::Edge; msg_labels=Dict{Interface, String}(), is_external_edge=false)
    # Generate DOT code for an edge
    dot = ""
    if edge.b != nothing
        b_id = findfirst(isequal(edge.b), edge.b.node.interfaces)
        b_label = "$b_id $(handle(edge.b))"
        b_object_id = objectid(edge.b.node)
    else
        b_object_id = "$(objectid(edge))2"
        dot *= "\t$b_object_id [shape=none, label=\"\", width=0.75, height=0.75]\n"
        b_label = ""
    end

    if edge.a != nothing
        a_id = findfirst(isequal(edge.a), edge.a.node.interfaces)
        a_label = "$a_id $(handle(edge.a))"
        a_object_id = objectid(edge.a.node)
    else
        a_object_id = "$(objectid(edge))1"
        dot *= "\t$a_object_id [shape=none, label=\"\", width=0.75, height=0.75]\n"
        a_label = ""
    end

    dot *= "\t$b_object_id -- $a_object_id "
    if is_external_edge
        dot *= "[taillabel=\"$(b_label)\", headlabel=\"$(a_label)\", style=\"dashed\" color=\"black\"]\n"
    else
        msg_label_a = haskey(msg_labels, edge.a) ? ""*msg_labels[edge.a]*"" : ""
        msg_label_b = haskey(msg_labels, edge.b) ? ""*msg_labels[edge.b]*"" : ""

        props = String[]
        if !isempty(b_label) || !isempty(msg_label_b)
            push!(props, "taillabel=<$(b_label) \n$(msg_label_b)>")
        end
        if !isempty(a_label) || !isempty(msg_label_a)
            push!(props, "headlabel=<$(a_label) \n$(msg_label_a)>")
        end
        push!(props, "label=<$(edge.variable.id)>")
        if !isempty(props)
            dot *= "[$(join(props,", "))]"
        end
    end
end

function dot2gif(dot_graph::AbstractString)
    # Generate SVG image from DOT graph
    validateGraphVizInstalled() # Show an error if GraphViz is not installed correctly
    proc = open(`dot -Tgif`, "r+")
    write(proc.in, dot_graph)
    close(proc.in)
    return read(proc.out, String)
end

function dot2svg(dot_graph::AbstractString)
    # Generate SVG image from DOT graph
    validateGraphVizInstalled() # Show an error if GraphViz is not installed correctly
    proc = open(`dot -Tsvg`, "r+")
    write(proc.in, dot_graph)
    close(proc.in)
    return read(proc.out, String)
end

function dot2png(dot_graph::AbstractString)
    # Generate PNG image from DOT graph
    validateGraphVizInstalled()
    proc = open(`dot -Tpng`, "r+")
    write(proc.in, dot_graph)
    close(proc.in)
    return read(proc.out, String)
end

function validateGraphVizInstalled()
    # Check if GraphViz is installed
    try
        (read(`dot -'?'`, String)[1:10] == "Usage: dot") || error()
    catch
        error("GraphViz is not installed correctly. Make sure GraphViz is installed. If you are on Windows, manually add the path to GraphViz to your path variable. You should be able to run 'dot' from the command line.")
    end
end

viewDotExternal(dot_graph::AbstractString) = (Sys.islinux() ? viewDotExternalInteractive(dot_graph::AbstractString) : viewDotExternalImage(dot_graph::AbstractString))

function viewDotExternalInteractive(dot_graph::AbstractString)
    # View a DOT graph in interactive viewer
    validateGraphVizInstalled() # Show an error if GraphViz is not installed correctly
    proc = open(`dot -Tx11`, "w")
    write(proc.in, dot_graph)
    close(proc.in)
end

function viewDotExternalImage(dot_graph::AbstractString)
    # Write the image to a file and open it with the default image viewer
    svg = dot2svg(dot_graph)
    filename = tempname()*".svg"
    open(filename, "w") do f
        write(f, svg)
    end
    viewFile(filename)
end

# Based on imgcat script provided by iTerm developers (working in iTerm v.2.9.x)
function viewDotIniTerm(dot_graph::AbstractString)
    png = dot2png(dot_graph)
    base64png = base64encode(png)
    sequence = ("$(Char(0x1b))]1337;File=size=$(length(png));inline=1:$base64png\a\n")
    print(sequence)
end


text:
# ---
# title: 993. Cousins in Binary Tree
# id: problem993
# author: Indigo
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Breadth-first Search
# link: <https://leetcode.com/problems/cousins-in-binary-tree/description/>
# hidden: true
# ---
#
# In a binary tree, the root node is at depth `0`, and children of each depth
# `k` node are at depth `k+1`.
#
# Two nodes of a binary tree are _cousins_ if they have the same depth, but have
# **different parents**.
#
# We are given the `root` of a binary tree with unique values, and the values
# `x` and `y` of two different nodes in the tree.
#
# Return `true` if and only if the nodes corresponding to the values `x` and `y`
# are cousins.
#
#
#
# **Example 1:
# ![](https://assets.leetcode.com/uploads/2019/02/12/q1248-01.png)**
#
#
#
#     Input: root = [1,2,3,4], x = 4, y = 3
#     Output: false
#
#
# **Example 2:
# ![](https://assets.leetcode.com/uploads/2019/02/12/q1248-02.png)**
#
#
#
#     Input: root = [1,2,3,null,4,null,5], x = 5, y = 4
#     Output: true
#
#
# **Example 3:**
#
# **![](https://assets.leetcode.com/uploads/2019/02/13/q1248-03.png)**
#
#
#
#     Input: root = [1,2,3,null,4], x = 2, y = 3
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree will be between `2` and `100`.
#   * Each node has a unique integer value from `1` to `100`.
#
#
## @lc code=start
using LeetCode

function is_cousins(root::TreeNode{Int}, x, y)
    q = NTuple{2, Union{Nothing, TreeNode}}[]
    push!(q, (root, nothing))
    while !isempty(q)
        n = length(q)
        rec_parent = TreeNode{Int}[]
        for i in 1:n
            cur, par = popfirst!(q)
            if cur.val == x || cur.val == y
                push!(rec_parent, par)
            end
            isnothing(cur.left) || push!(q, (cur.left, cur))
            isnothing(cur.right) || push!(q, (cur.right, cur))
        end
        length(rec_parent) == 0 && continue
        length(rec_parent) == 1 && return false
        length(rec_parent) == 2 && return rec_parent[1] != rec_parent[2]
    end
end
## @lc code=end

ext: jl


text:
# Someb of this is a translation of Canary.jl's metric tests written by Jeremy
# Kozdon.

@testset "metrics" begin
    TAs = ((Float64,  Array),
           (Float32,  Array),
           # (BigFloat, Array)
          )
    if CUDA.has_cuda_gpu()
        TAs = (TAs..., (Float32, CuArray))
    end

    for (T, A) in TAs

        @testset "1D" begin
            cell = LobattoCell{T, A}(4);
            grid = brickgrid(cell, ([-1, 0, 10],))
            g, J = adapt(Array, components(metrics(grid)))
            g, = components(g)

            @test all(g[:, 1] .≈ 2)
            @test all(g[:, 2] .≈ 1//5)

            @test all(J[:, 1] .≈ 1//2)
            @test all(J[:, 2] .≈ 5)

            n, sJ = adapt(Array, components(facemetrics(grid)))
            n, = components(n)

            @test all(n[1, :] .≈ -1)
            @test all(n[2, :] .≈ 1)
            @test all(sJ[:] .≈ 1)
        end

        @testset "2D" begin
            f(x) = SA[9 * x[1] - (1 + x[1]) * x[2]^2 +
                      (x[1] - 1)^2 * (1 - x[2]^2 + x[2]^3),
                      10 * x[2] + x[1] * x[1]^3 * (1 - x[2]) +
                      x[1]^2 * x[2] * (1 + x[2])]
            f₁₁(x) = 7 + x[2]^2 - 2 * x[2]^3 + 2 * x[1] * (1 - x[2]^2 + x[2]^3)
            f₁₂(x) = -2 * (1 + x[1]) * x[2] +
                     (-1 + x[1])^2 * x[2] * (-2 + 3 * x[2])
            f₂₁(x) = -4 * x[1]^3 * (-1 + x[2]) + 2 * x[1] * x[2] * (1 + x[2])
            f₂₂(x) = 10 - x[1]*x[1]^3 + x[1]^2 * (1 + 2 * x[2])

            L = 10
            M = 12
            cell = LobattoCell{T, A}(L, M)
            Δx₁ = 1//5
            Δx₂ = 1//4

            unwarpedgrid = brickgrid(cell, (-1:Δx₁:1, -1:Δx₂:1))
            x̂ = points(unwarpedgrid)
            ĥ = fieldarray(SMatrix{2, 2, T, 4},
                           (f₁₁.(x̂).*(Δx₁/2), f₂₁.(x̂).*(Δx₁/2),
                            f₁₂.(x̂).*(Δx₂/2), f₂₂.(x̂).*(Δx₂/2)))
            Ĵ = adapt(Array, det.(ĥ))
            ĝ = adapt(Array, inv.(ĥ))

            grid = brickgrid(f, cell, (-1:Δx₁:1, -1:Δx₂:1))
            g, J = map(x->adapt(Array, x), components(metrics(grid)))

            @test J ≈ Ĵ
            @test all(g .≈ ĝ)

            n, sJ = map(x->adapt(Array, x), components(facemetrics(grid)))

            @test all(norm.(n) .≈ 1)

            a = n .* sJ
            a₁, a₂ = components(a)
            ĥ₁₁, ĥ₂₁, ĥ₁₂, ĥ₂₂ = map(x->adapt(Array, x),
                                     components(reshape(ĥ, size(cell)..., :)))

            @test a₁[1:M, :] ≈ -ĥ₂₂[1, :, :]
            @test a₂[1:M, :] ≈  ĥ₁₂[1, :, :]

            @test a₁[M+1:2M, :] ≈  ĥ₂₂[end, :, :]
            @test a₂[M+1:2M, :] ≈ -ĥ₁₂[end, :, :]

            @test a₁[2M+1:2M+L, :] ≈  ĥ₂₁[:, 1, :]
            @test a₂[2M+1:2M+L, :] ≈ -ĥ₁₁[:, 1, :]

            @test a₁[2M+L+1:2M+2L, :] ≈ -ĥ₂₁[:, end, :]
            @test a₂[2M+L+1:2M+2L, :] ≈  ĥ₁₁[:, end, :]
        end

        @testset "2D Constant Preserving" begin
            f(x) = SA[9 * x[1] - (1 + x[1]) * x[2]^2 +
                      (x[1] - 1)^2 * (1 - x[2]^2 + x[2]^3),
                      10 * x[2] + x[1] * x[1]^3 * (1 - x[2]) +
                      x[1]^2 * x[2] * (1 + x[2])]

            L = 3
            M = 4
            cell = LobattoCell{T, A}(L, M)
            Δx₁ = 1//1
            Δx₂ = 1//2

            grid = brickgrid(f, cell, (-1:Δx₁:1, -1:Δx₂:1))
            g, J = map(x->adapt(Array, x), components(metrics(grid)))
            D₁, D₂ = adapt(Array, derivatives(cell))

            @test norm((D₁*(J .* g.:1) + D₂*(J .* g.:2)), Inf) < 100*eps(T)
            @test norm((D₁*(J .* g.:3) + D₂*(J .* g.:4)), Inf) < 100*eps(T)
        end

        @testset "3D" begin
            f(x) = SA[3x[1] + x[2]/5 + x[3]/10 + x[1]*x[2]^2*x[3]^3/3,
                      4x[2] + x[1]^3*x[2]^2*x[3]/4,
                      2x[3] + x[1]^2*x[2]*x[3]^3/2]
            f₁₁(x) = 3*oneunit(eltype(x)) + x[2]^2*x[3]^3/3
            f₁₂(x) = oneunit(eltype(x))/5 + 2*x[1]*x[2]*x[3]^3/3
            f₁₃(x) = oneunit(eltype(x))/10 + 3*x[1]*x[2]^2*x[3]^2/3
            f₂₁(x) = 3*x[1]^2*x[2]^2*x[3]/4
            f₂₂(x) = 4*oneunit(eltype(x)) + 2*x[1]^3*x[2]*x[3]/4
            f₂₃(x) = x[1]^3*x[2]^2/4
            f₃₁(x) = 2*x[1]*x[2]*x[3]^3/2
            f₃₂(x) = x[1]^2*x[3]^3/2
            f₃₃(x) = 2*oneunit(eltype(x)) + 3*x[1]^2*x[2]*x[3]^2/2

            L = 6
            M = 8
            N = 7
            cell = LobattoCell{T, A}(L, M, N)
            Δx₁ = 1//1
            Δx₂ = 2//1
            Δx₃ = 2//3

            unwarpedgrid = brickgrid(cell, (-1:Δx₁:1, -1:Δx₂:1, -1:Δx₃:1))
            x̂ = points(unwarpedgrid)
            ĥ = fieldarray(SMatrix{3, 3, T, 9},
                           (
                            f₁₁.(x̂).*(Δx₁/2),
                            f₂₁.(x̂).*(Δx₁/2),
                            f₃₁.(x̂).*(Δx₁/2),
                            f₁₂.(x̂).*(Δx₂/2),
                            f₂₂.(x̂).*(Δx₂/2),
                            f₃₂.(x̂).*(Δx₂/2),
                            f₁₃.(x̂).*(Δx₃/2),
                            f₂₃.(x̂).*(Δx₃/2),
                            f₃₃.(x̂).*(Δx₃/2),
                           ))
            Ĵ = adapt(Array, det.(ĥ))
            ĝ = adapt(Array, inv.(ĥ))

            grid = brickgrid(f, cell, (-1:Δx₁:1, -1:Δx₂:1, -1:Δx₃:1))
            g, J = map(x->adapt(Array, x), components(metrics(grid)))

            @test J ≈ Ĵ
            @test all(g .≈ ĝ)

            n, sJ = map(x->adapt(Array, x), components(facemetrics(grid)))

            @test all(norm.(n) .≈ 1)

            a = n .* sJ
            a₁, a₂, a₃ = components(a)
            ĝ₁₁, ĝ₂₁, ĝ₃₁, ĝ₁₂, ĝ₂₂, ĝ₃₂, ĝ₁₃, ĝ₂₃, ĝ₃₃ =
                map(x->adapt(Array, x),
                    components(reshape(ĝ, size(cell)..., :)))
            Ĵ = reshape(Ĵ, size(cell)..., :)

            @test reshape(a₁[      1:M*N, :], M, N, :) ≈ -Ĵ[  1, :, :, :] .* ĝ₁₁[  1, :, :, :]
            @test reshape(a₁[M*N+1:2*M*N, :], M, N, :) ≈  Ĵ[end, :, :, :] .* ĝ₁₁[end, :, :, :]
            @test reshape(a₁[ 2*M*N        .+ (1:L*N), :], L, N, :) ≈ -Ĵ[:,   1, :, :] .* ĝ₂₁[:,   1, :, :]
            @test reshape(a₁[(2*M*N + L*N) .+ (1:L*N), :], L, N, :) ≈  Ĵ[:, end, :, :] .* ĝ₂₁[:, end, :, :]
            @test reshape(a₁[(2*M*N + 2*L*N      ) .+ (1:L*M), :], L, M, :) ≈ -Ĵ[:, :,   1, :] .* ĝ₃₁[:, :,   1, :]
            @test reshape(a₁[(2*M*N + 2*L*N + L*M) .+ (1:L*M), :], L, M, :) ≈  Ĵ[:, :, end, :] .* ĝ₃₁[:, :, end, :]

            @test reshape(a₂[      1:M*N, :], M, N, :) ≈ -Ĵ[  1, :, :, :] .* ĝ₁₂[  1, :, :, :]
            @test reshape(a₂[M*N+1:2*M*N, :], M, N, :) ≈  Ĵ[end, :, :, :] .* ĝ₁₂[end, :, :, :]
            @test reshape(a₂[ 2*M*N        .+ (1:L*N), :], L, N, :) ≈ -Ĵ[:,   1, :, :] .* ĝ₂₂[:,   1, :, :]
            @test reshape(a₂[(2*M*N + L*N) .+ (1:L*N), :], L, N, :) ≈  Ĵ[:, end, :, :] .* ĝ₂₂[:, end, :, :]
            @test reshape(a₂[(2*M*N + 2*L*N      ) .+ (1:L*M), :], L, M, :) ≈ -Ĵ[:, :,   1, :] .* ĝ₃₂[:, :,   1, :]
            @test reshape(a₂[(2*M*N + 2*L*N + L*M) .+ (1:L*M), :], L, M, :) ≈  Ĵ[:, :, end, :] .* ĝ₃₂[:, :, end, :]

            @test reshape(a₃[      1:M*N, :], M, N, :) ≈ -Ĵ[  1, :, :, :] .* ĝ₁₃[  1, :, :, :]
            @test reshape(a₃[M*N+1:2*M*N, :], M, N, :) ≈  Ĵ[end, :, :, :] .* ĝ₁₃[end, :, :, :]
            @test reshape(a₃[ 2*M*N        .+ (1:L*N), :], L, N, :) ≈ -Ĵ[:,   1, :, :] .* ĝ₂₃[:,   1, :, :]
            @test reshape(a₃[(2*M*N + L*N) .+ (1:L*N), :], L, N, :) ≈  Ĵ[:, end, :, :] .* ĝ₂₃[:, end, :, :]
            @test reshape(a₃[(2*M*N + 2*L*N      ) .+ (1:L*M), :], L, M, :) ≈ -Ĵ[:, :,   1, :] .* ĝ₃₃[:, :,   1, :]
            @test reshape(a₃[(2*M*N + 2*L*N + L*M) .+ (1:L*M), :], L, M, :) ≈  Ĵ[:, :, end, :] .* ĝ₃₃[:, :, end, :]
        end

        @testset "3D Constant Preserving" begin
            f(x) = SA[3x[1] + x[2]/5 + x[3]/10 + x[1]*x[2]^2*x[3]^3/3,
                      4x[2] + x[1]^3*x[2]^2*x[3]/4,
                      2x[3] + x[1]^2*x[2]*x[3]^3/2]

            L = 3
            M = 4
            N = 2
            cell = LobattoCell{T, A}(L, M, N)
            Δx₁ = 1//1
            Δx₂ = 1//2
            Δx₃ = 2//3

            grid = brickgrid(f, cell, (-1:Δx₁:1, -1:Δx₂:1, -1:Δx₃:1))
            g, J = map(x->adapt(Array, x), components(metrics(grid)))
            D₁, D₂, D₃ = adapt(Array, derivatives(cell))

            @test norm((D₁*(J .* g.:1) + D₂*(J .* g.:2) + D₃*(J .* g.:3)),
                       Inf) < 100*eps(T)
            @test norm((D₁*(J .* g.:4) + D₂*(J .* g.:5) + D₃*(J .* g.:6)),
                       Inf) < 100*eps(T)
            @test norm((D₁*(J .* g.:7) + D₂*(J .* g.:8) + D₃*(J .* g.:9)),
                       Inf) < 100*eps(T)
        end
    end
end


text:
"""
    simulate(
        tini, tfin, dss, controls, aircraft, atmosphere, gravity;
        solver=TSit5(), solve_args=Dict()
    )

"""
function simulate(tini, tfin, dss, args...; kwargs...)

    sol = _simulate(tini, tfin, dss, args...; solver, solve_args)

    df = DataFrame(sol')
    rename!(df, get_x_names(dss))
    df[!, :time] = sol.t

    return df
end
function _simulate(tini, tfin, dss, controls, aircraft, atmosphere, gravity;
    solver=Tsit5(), solve_args=Dict()
    )

    tspan = (tini, tfin)
    p = [dss, controls, aircraft, atmosphere, gravity]

    x0 = get_x(dss)
    prob = ODEProblem{false}(f, x0, tspan, p)
    return solve(prob, solver; solve_args...)
end

function f(x, p, t)
    dss = p[1]
    controls = p[2]
    aircraft = p[3]
    atmosphere = p[4]
    gravity = p[5]

    controls_arr = get_value.(controls, t)

    dss = typeof(dss)(x)

    atmosphere = atmosphere(get_height(dss))

    dssd, outputs = f(time, dss, controls_arr, aircraft, atmosphere, gravity)

    return get_xdot(dssd)
end

function f(
    time,
    dss::DSState,
    controls,
    aircraft::Aircraft,
    atmosphere::Atmosphere,
    gravity::Gravity
    )

    ac = aircraft

    mass = get_mass(ac)
    xcg = get_xcg_mac(ac)
    inertia = get_inertia_tensor(ac)

    S = get_surface(ac)
    c = get_chord(ac)
    b = get_wing_span(ac)

    vt = get_tas(dss)
    α = get_α(dss) * RAD2DEG
    β = get_β(dss) * RAD2DEG
    ψ, θ, ϕ = get_euler_angles(dss)
    p, q, r = get_ang_vel_body(dss)
    height = get_height(dss)
    pow = get_engine_power(dss)

    T = get_temperature(atmosphere)
    ρ = get_density(atmosphere)
    a = get_sound_velocity(atmosphere)
    p = get_pressure(atmosphere)

    # Air data computer
    # TODO: this will be inside AeroState
    amach = vt / a  # mach number
    qbar = 0.5 * ρ * vt^2  # dynamic pressure

    # Calculate forces and moments
    # Propulsion
    Tx, Ty, Tz, LT, MT, NT = calculate_prop_forces_moments(ac, dss, amach, controls)
    h = calculate_prop_gyro_effects(ac)

    # Engine dynamic model
    thtl = controls[1]
    pdot = calculate_pdot(ac, thtl, pow)

    # Aerodynamics
    Fax, Fay, Faz, La, Ma, Na = calculate_aero_forces_moments(ac, dss, controls, xcg, qbar, S, b, c)
    # Gravity
    Fgx, Fgy, Fgz = get_gravity_body(gravity, θ, ϕ) .* mass

    # Total forces & moments
    Fx = Fgx + Fax + Tx
    Fy = Fgy + Fay + Ty
    Fz = Fgz + Faz + Tz

    L = La
    M = Ma
    N = Na

    forces = [Fx, Fy, Fz]
    moments = [L, M, N]

    dssd = state_eqs(
        dss, time, mass, inertia, forces, moments, h, pdot
    )

    # Outputs
    gravity_down = get_gravity_accel(gravity)
    outputs = calculate_outputs(dss, amach, qbar, S, mass, gravity_down, [Fax, Fay, Faz], [Tx, Ty, Tz])

    return dssd, outputs

end

function calculate_outputs(dss::DSState, amach, qbar, S, mass, g, Fa, Fp)

    outputs = Array{Float64}(undef, 7)

    α = get_α(dss) * RAD2DEG
    p, q, r = get_ang_vel_body(dss)

    Fax, Fay, Faz = Fa
    Tx, Ty, Tz = Fp

    qbarS = qbar * S

    rmqs = qbarS / mass

    ax = (Fax + Tx) / mass
    ay = (Fay + Ty) / mass
    az = (Faz + Tz) / mass

    along = ax / g
    an = -az / g
    alat = ay / g

    outputs[1] = an
    outputs[2] = alat
    outputs[3] = along
    outputs[4] = qbar
    outputs[5] = amach
    outputs[6] = q
    outputs[7] = α

    return outputs
end


text:
"""
Generate structures using the `generate_structs()` function from InfrastructureSystems.
"""
# SRC_DIR = abspath(dirname(Base.find_package("SLiDE")))
SRC_DIR = joinpath(SLIDE_DIR, "src")

descriptor_file = joinpath(SRC_DIR, "descriptors", "check_data_structs.json")
output_dir = joinpath(SRC_DIR, "parse", "generated_check")

IS.generate_structs(descriptor_file, output_dir; print_results = true)

descriptor_file = joinpath(SRC_DIR, "descriptors", "edit_data_structs.json")
output_dir = joinpath(SRC_DIR, "parse", "generated_edit")

IS.generate_structs(descriptor_file, output_dir; print_results = true)

descriptor_file = joinpath(SRC_DIR, "descriptors", "load_data_structs.json")
output_dir = joinpath(SRC_DIR, "parse", "generated_load")

IS.generate_structs(descriptor_file, output_dir; print_results = true)

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
module FastbackTests

using Test
using Dates
using Fastback

include("permute.jl")
include("bidask.jl")
include("position.jl")
include("collectors.jl")
include("account.jl")

end


text:
# Definition of CuckooFilter and associated methods

"""
    CuckooFilter{F}

A cuckoo filter. `F` is the type of fingerprint used by the filter (defaults
to `UInt8`).

# Examples

```jldoctest; setup = :(using CuckooFilters)
julia> filter = CuckooFilter{AbstractVector{UInt8}}();

julia> b"hello, world!" ∈ filter
false

julia> insert!(filter, b"hello, world!")
true

julia> b"hello, world!" ∈ filter
true

julia> delete!(filter, b"hello, world!")
true

julia> b"hello, world!" ∈ filter
false
```
"""
struct CuckooFilter{T,F}
    buckets :: Array{Bucket{F},1}
    max_bucket_kicks :: Int64
    max_bucket_size :: Int64
end

CuckooFilter{T}() where T = CuckooFilter{T, fingerprint_type(T)}()

function CuckooFilter{T,F}(;
    n_buckets = 1000,
    max_bucket_kicks = 100,
    max_bucket_size = 4
) where {T,F}
    buckets = [Bucket{F}() for i = 1:n_buckets]
    CuckooFilter{T,F}(buckets, 100, 4)
end

function getbucket(filter, i)
    n_buckets = length(filter.buckets)
    i = (i % n_buckets) + 1
    return filter.buckets[i]
end

insert_unique!(filter::CuckooFilter, key) = insert!(filter, key, true)

function Base.insert!(
    filter::CuckooFilter{T},
    key::T,
    unique = false
) where T
    if unique && key ∈ filter
        return false
    end

    f = fingerprint(key)
    h₁ = hash(key)

    # Probe first bucket to see if it has an empty entry
    b₁ = getbucket(filter, h₁)
    if !isfull(b₁, filter.max_bucket_size)
        insert!(b₁, f)
        return true
    end

    # Probe second bucket to see if it has an empty entry
    h₂ = h₁ ⊻ hash(f)
    b₂ = getbucket(filter, h₂)
    if !isfull(b₂, filter.max_bucket_size)
        insert!(b₂, f)
        return true
    end

    # If both buckets are full, we need to relocate an existing
    # item from one of the two buckets and move it to another
    # bucket
    b, h = rand(((h₁, b₁), (h₂, b₂)))
    for i = 1:filter.max_bucket_kicks
        original_f = f

        # Swap the fingerprint f with a random fingerprint from
        # the bucket
        rand_idx = rand(1:length(b))
        f = b[rand_idx]
        b[rand_idx] = original_f

        # Attempt to insert the fingerprint into a new bucket
        h = h ⊻ hash(f)
        b = getbucket(f, h)
        if !isfull(b)
            insert!(b, f)
            return true
        end
    end

    # If we've reached this point, we've gone through the maximum
    # number of possible iterations to try and kick fingerprints
    # through the cuckoo filter. We treat the hash table as full.
    error("Cuckoo filter is full")
end

function Base.delete!(filter::CuckooFilter{T}, key::T) where T
    # Delete the key from its first possible bucket if it is inside
    # of that bucket
    f = fingerprint(key)
    h₁ = hash(key)
    b₁ = getbucket(filter, h₁)
    if delete!(b₁, f)
        return true
    end

    # If the first deletion failed, then we attempt to delete the
    # key from its second bucket
    h₂ = h₁ ⊻ hash(f)
    b₂ = getbucket(filter, h₂)
    delete!(b₂, f)
end

function Base.in(key::T, filter::CuckooFilter{T}) where T
    # Check whether the key's fingerprint is in its first possible
    # bucket
    f = fingerprint(key)
    h₁ = hash(key)
    if f ∈ getbucket(filter, h₁)
        return true
    end

    # Check whether the fingerprint is in its second possible bucket
    h₂ = h₁ ⊻ hash(f)
    f ∈ getbucket(filter, h₂)
end


text:
struct QueryableFilter <: Queryable
    source
    filter_func
    filter_expr
    getiterator
end

function QueryOperators.filter(source::Queryable, filter_func::Function, filter_expr::Expr)
    return QueryableFilter(source, filter_func, filter_expr, source.getiterator)
end


text:
using AutomotivePOMDPs
using Test

include("collisionchecker_tests.jl")
include("crosswalk_tests.jl")
include("urbanpomdp_tests.jl")

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
module BenchLoadBalancing

using BenchmarkTools
using Folds
using FoldsThreads

const SUITE = BenchmarkGroup()
const DATA_LENGTH = 2^10

function unbalanced_work(nworks, data_length)
    nworks == 0 && return nowork(_) = 0
    space = cld(data_length, nworks)
    function work(i)
        n = 0
        if mod(i, space) == 0
            # Sping for 100 μs
            t = time_ns() + 100_0000
            while t > time_ns()
                n += 1
            end
        end
        return n
    end
end

for nworks in [0, 5, 10, 30, 35]
    f = unbalanced_work(nworks, DATA_LENGTH)
    xs = 1:DATA_LENGTH
    SUITE[(nworks = nworks,)] =
        @benchmarkable(Folds.sum($f, $xs, WorkStealingEx(basesize = 1)))
end

end  # module
BenchLoadBalancing.SUITE


text:
function solve(d, n=1)
    l = collect(0:255)
    i = 0
    skip = 0
    for _ in 1:n
        for j in d
            l = circshift(l, -i)
            reverse!(l, 1, j)
            l = circshift(l, i)
            i += j + skip
            skip += 1
        end
    end
    return l
end

# Part 1
d = readdlm("input.txt", ',', Int)
l = solve(d)
println(l[1] * l[2])

# Part 2
d = vcat(map(Int, read("input.txt")[1:end-1]), [17,31,73,47,23])
l = solve(d, 64)
l = mapslices(x -> reduce(⊻, x), reshape(l, (16,16)), 1)
println(mapreduce(hex, *, l))


text:
function loglike(s::State{A1, A2, A3}, y::Vector{MA},
                 m::Vector{BitArray{2}}, c::Constants) where {A1, A2, A3, MA}

  sig = sqrt.(s.sig2)
  noisy_sd = sqrt(c.noisy_var)

  # reshape mu[z] to have size (1 x 1 x L[z])
  mu0 = reshape(-cumsum(s.delta0), 1, 1, c.L[0])
  mu1 = reshape(cumsum(s.delta1), 1, 1, c.L[1])

  # Z: J x K
  # H: J x K
  # v: K
  Z = compute_Z(s.v, s.H, tau=c.tau, use_stickbreak=c.use_stickbreak)
  Z_rs = reshape(Z, 1, c.J, c.K)

  # ll = zero(s.alpha)
  ll = 0.0
  for i in 1:c.I
    mi = m[i]
    Ni = size(y[i], 1)

    # Ni x J x Lz
    yi = reshape(y[i], Ni, c.J, 1)

    # Ni x J x 1
    logmix_L0 = ADVI.lpdf_gmm(yi, mu0, sig[i], s.eta0[i:i, :, :],
                              dims=3, dropdim=false)
    logmix_L1 = ADVI.lpdf_gmm(yi, mu1, sig[i], s.eta1[i:i, :, :],
                              dims=3, dropdim=false)
    @assert size(logmix_L1) == (Ni, c.J, 1) == size(logmix_L1)

    # Ni x J x K -> Ni x K
    Z_mix = ADVI.sumdd(Z_rs .* logmix_L1 + (1 .- Z_rs) .* logmix_L0, dims=2)

    # Ni x K
    f = Z_mix .+ log.(s.W[i:i, :])

    # Ni - dimensional
    lli_pre = ADVI.logsumexpdd(f, dims=2)

    # mix with noisy
    lli_quiet = lli_pre .+ log1p(-s.eps[i])
    lli_noisy = ADVI.sumdd(ADVI.lpdf_normal.(y[i], 0., noisy_sd), dims=2) .+ log(s.eps[i])

    # Ni - dimensional
    lli = ADVI.logsumexpdd(ADVI.stack(lli_quiet, lli_noisy), dims=-1)

    # p(m | y)
    pm_i = prob_miss(y[i][mi], c.beta[i]...)
    logprob_mi_given_yi = sum(log.(pm_i))

    # add to ll
    fac = c.N[i] / Ni
    ll += (sum(lli) + logprob_mi_given_yi) * fac
  end

  return ll
end

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using Documenter, Plots, PlotDocs

for be in (:gr, :pyplot, :plotlyjs, :pgfplots)
    generate_markdown(be)
end
save_attr_html_files()
create_support_tables()

makedocs(modules=[Plots], doctest=false)

# makedocs(
#     format = :html,
#     sitename = "Plots",
#     authors = "Thomas Breloff",
#     pages = Any[
#         "Home" => "index.md",
#         "Getting Started" => [
#             "Installation" => "install.md",
#             "Basics" => "basics.md",
#             "Tutorial" => "tutorial.md",
#         ],
#         "Manual" => [
#             "Input Data" => "input_data.md",
#             "Output" => "output.md",
#             "Attributes" => "attributes.md",
#             "Layouts" => "layouts.md",
#             "Recipes" => "recipes.md",
#             "Colors" => "colors.md",
#             "Animations" => "animations.md",
#             "Backends" => "backends.md",
#             "Supported Attributes" => "supported.md",
#         ],
#         "Learning" => "learning.md",
#         "Contributing" => "contributing.md",
#         "Ecosystem" => "ecosystem.md",
#         "Advanced Topics" => [
#             "Internals" => "pipeline.md",
#         ],
#         "Examples" => [
#             "GR" => "examples/gr.md",
#             "PyPlot" => "examples/pyplot.md",
#             "PlotlyJS" => "examples/plotlyjs.md",
#             "PGFPlots" => "examples/pgfplots.md",
#             "UnicodePlots" => "examples/unicodeplots.md",
#             "InspectDR" => "examples/inspectdr.md",
#         ],
#     ],
#     )

# deploydocs(
#     repo = "github.com/JuliaPlots/PlotDocs.jl.git",
#     target = "build",
#     deps = nothing,
#     make = nothing,
# )

deploydocs(
    deps = Deps.pip("mkdocs==0.17.5",
        "mkdocs-material==2.9.4" ,"python-markdown-math", "pygments", "pymdown-extensions"),
    repo   = "github.com/JuliaPlots/PlotDocs.jl.git",
    julia  = "1.0",
    osname = "linux"
)


text:
using stoosc
rc("font", family="arial")

function params_estimate_arnold_tongue!(params)

    # Extract function parameters
    trajectories = params["trajectories"]
    parameter_names = params["parameter_names"]
    parameter_values = params["parameter_values"]
    n_samples = params["n_samples"]
    model_name = params["model_name"]
    I_max = params["I_max"]
    problem_type = params["problem_type"]
    i_variable = params["i_variable"]
    period_range = params["period_range"]
    filename = params["filename"]

    # Create a model
    model = load_model(model_name, problem_type)
    model = set_solver(model, saveat=0.01, maxiters=1e10)
    model = set_timespan(model, (0.0, 200.0))
    outfun = sol -> hcat(sol[i_variable, :])
    model = set_output(model, outfun)
    model = set_parameter_value(model, parameter_names, parameter_values)

    # Set periods and amplitudes to scan
    periods = range(period_range[1], period_range[2], length=n_samples)
    # periods = range(0.8, 1.2, length=n_samples)
    amplitudes = [0.0, I_max]

    # Estimate arnold tongue
    arnold_tongue = estimate_binary_arnold_tongue(model, periods, amplitudes;
        trajectories=trajectories, i_solution=1, n_lags=1000,
        input_parameter_name="I", min_time=nothing, max_time=nothing,
        show_progress=true, period_error=0.02, boundary_error=I_max/n_samples,
        maximal_error_attempts=5)

    # Save the estimated arnold tongue into the params dictionary
    params["arnold_tongue"] = arnold_tongue

    # Save the estimated arnold tongue into a file
    if !isnothing(filename)
        save_data(arnold_tongue, filename)
    end

end

params = Dict(
    "trajectories" => 10,  # number of trajectories to calculate population mean
    "parameter_names" => ["σ"],  # parameter names
    "parameter_values" => [0.005],  # parameter values
    "n_samples" => 100,  # number of samples in the period range
    "period_range" => [0.75, 1.25],  # period range (min, max)
    "model_name" => "model01",  # model name (model01, model02)
    "I_max" => 0.3,  # maximal input amplitude to estimate
    "problem_type" => "sde",  # "ode", "sde", or "jump"
    "i_variable" => 3,  # varialbe to use to estimate entrainment
    "filename" => "test_arnold_tongue.csv",  # save the result in this file
)
params_estimate_arnold_tongue!(params)

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
immutable Uniform <: ContinuousUnivariateDistribution
    a::Float64
    b::Float64

    function Uniform(a::Real, b::Real)
	   a < b || error("Uniform: a must be less than b")
	   new(float64(a), float64(b))
    end

    Uniform() = new(0.0, 1.0)
end

@distr_support Uniform d.a d.b

#### Parameters

params(d::Uniform) = (d.a, d.b)

location(d::Uniform) = d.a
scale(d::Uniform) = d.b - d.a

#### Statistics

mean(d::Uniform) = middle(d.a, d.b)
median(d::Uniform) = mean(d)
mode(d::Uniform) = mean(d)
modes(d::Uniform) = Float64[]

var(d::Uniform) = (w = d.b - d.a; w^2 / 12.0)

skewness(d::Uniform) = 0.0
kurtosis(d::Uniform) = -1.2

entropy(d::Uniform) = log(d.b - d.a)

#### Evaluation

pdf(d::Uniform, x::Float64) = insupport(d, x) ? 1.0 / (d.b - d.a) : 0.0
logpdf(d::Uniform, x::Float64) = insupport(d, x) ? -log(d.b - d.a) : -Inf

function cdf(d::Uniform, x::Float64)
    (a, b) = params(d)
    x <= a ? 0.0 :
    x >= d.b ? 1.0 : (x - a) / (b - a)
end

function ccdf(d::Uniform, x::Float64)
    (a, b) = params(d)
    x <= a ? 1.0 :
    x >= d.b ? 0.0 : (b - x) / (b - a)
end

quantile(d::Uniform, p::Float64) = d.a + p * (d.b - d.a)
cquantile(d::Uniform, p::Float64) = d.b + p * (d.a - d.b)

function mgf(d::Uniform, t::Real)
    (a, b) = params(d)
    u = 0.5 * (b - a) * t
    u == zero(u) && return one(u)
    v = 0.5 * (a + b) *t
    exp(v) * (sinh(u) / u)
end

function cf(d::Uniform, t::Real)
    (a, b) = params(d)
    u = 0.5 * (b - a) * t
    u == zero(u) && return complex(one(u))
    v = 0.5 * (a + b) * t
    exp(im * v) * (sin(u) / u)
end

#### Evaluation

rand(d::Uniform) = d.a + (d.b - d.a) * rand()

#### Fitting

function fit_mle{T <: Real}(::Type{Uniform}, x::Vector{T})
    if isempty(x)
        throw(ArgumentError("x cannot be empty."))
    end

    xmin = xmax = x[1]
    for i = 2:length(x)
        xi = x[i]
        if xi < xmin
            xmin = xi
        elseif xi > xmax
            xmax = xi
        end
    end

    Uniform(xmin, xmax)
end


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

# %% [markdown]
# # ベイズハッキング
#
# * 黒木玄
# * 2021-12-04
# * [nbviewer](https://nbviewer.org/github/genkuroki/public/blob/main/0025/Bayes%20hacking.ipynb)
#
# 関連ノートブック
#
# * [止め方で結果が変わる？](https://nbviewer.org/github/genkuroki/public/blob/main/0025/%E6%AD%A2%E3%82%81%E6%96%B9%E3%81%A7%E7%B5%90%E6%9E%9C%E3%81%8C%E5%A4%89%E3%82%8F%E3%82%8B%EF%BC%9F.ipynb)
# * [Clopper-Pearsonの信頼区間](https://nbviewer.org/github/genkuroki/public/blob/main/0025/Clopper-Pearson%20confidence%20interval.ipynb)

# %%
using Distributions
using StatsPlots
using StatsFuns: logit, logistic
using Roots
using StatsBase: ecdf
using Memoization
using SpecialFunctions

using Base64
showimg(mime, fn; tag="img") = open(fn) do f
    base64 = base64encode(f)
    display("text/html", """<$tag src="data:$mime;base64,$base64" />""")
end

safemul(x, y) = iszero(x) ? x : x*y
x ⪅ y = x < y || x ≈ y

# %%
showimg("image/jpeg", "IMG_2864.jpg"; tag="img width=50%")

# %%
showimg("image/jpeg", "IMG_2865.jpg"; tag="img width=60%")

# %%
"""
事前分布 Beta(a, b) のもとでn回中ｋ回成功の場合の事後分布は Beta(k+a, n-k+b) になる。
デフォルトでは a = b = 1 の一様事前分布を採用することにする。
"""
posterior_dist(n, k; a=1, b=a) = Beta(k + a, n - k + b)

"""pvalue_bayes(n, p, k; a, b) < α ⇔ p ∉ (ベイズ版100(1-α)%信用区間)"""
@memoize function pvalue_bayes(n, p, k; a = 1, b = a)
    posterior = posterior_dist(n, k; a, b)
    v0 = logpdf(posterior, p)
    f(t) = logpdf(posterior, logistic(t)) - v0
    m = params(posterior) |> ((α, β),) -> (α - 1)/(α + β - 2)
    if m ≤ 0
        s = ccdf(posterior, p)
    elseif m ≥ 1
        s = cdf(posterior, p)
    elseif p ≤ m
        q = logistic(find_zero(f, min(50, logit(m) + 1)))
        s = cdf(posterior, p) + ccdf(posterior, q)
    else
        q = logistic(find_zero(f, max(-50, logit(m) - 1)))
        s = cdf(posterior, q) + ccdf(posterior, p)
    end
    min(1, s)
end

"""
Bernoulli(p)分布の独立試行をNまで行うとき、
pがベイズ版信用区間から外れるまでの試行回数をL個計算
"""
function bayeshacking(N, p; α = 0.05, L = 10^4, a = 1, b = a)
    numtrials = fill(N + 1, L) # N＋１で埋め尽くされた長さLのベクトル
    Threads.@threads for i in 1:L # 並列のでL回のシミュレーションを実行
        k = 0
        for n in 1:N # ベルヌイ試行回数nの最大値はN
            k += rand(Bernoulli(p)) # kに1または0を足す。1を足す確率はp.
            if pvalue_bayes(n, p, k; a, b) < α # もしもpの値がベイズ信用区間を外れたら
                numtrials[i] = n  # n回で目的を達成したことを記録に残す
                break
            end
        end
    end
    numtrials
end

# %%
numtrials_bhack = bayeshacking(10^3, 0.5)
[(n, ecdf(numtrials_bhack)(n)) for n in (10, 100, 1000)]

# %%
plot(n -> ecdf(numtrials_bhack)(n), 0, 10^3; label="", xtick=0:100:1000)
title!("ecdf(number of trials for Bayes hacking)"; titlefontsize=12)

# %%
"""通常のP値函数(exact版)"""
@memoize function pvalue_exact(n, p, k)
    bin = Binomial(n, p)
    p0 = pdf(bin, k)
    s = sum(pdf(bin, j) for j in support(bin) if pdf(bin, j) ⪅ p0)
    min(1, s)
end

"""通常のP値函数(正規分布近似)"""
@memoize function pvalue_normal(n, p, k)
    bin = Binomial(n, p)
    normal = Normal(mean(bin), std(bin))
    min(1, 2cdf(normal, k), 2ccdf(normal, k))
end

"""
Bernoulli(p)分布の独立試行をNまで行うとき、
P値がα未満になるまでの試行回数をL個計算
"""
function phacking(pvalue_func, N, p; α = 0.05, L = 10^4)
    numtrials = fill(N + 1, L) # N＋１で埋め尽くされた長さLのベクトル
    Threads.@threads for i in 1:L # 並列のでL回のシミュレーションを実行
        k = 0
        for n in 1:N # ベルヌイ試行回数nの最大値はN
            k += rand(Bernoulli(p)) # kに1または0を足す。1を足す確率はp.
            if pvalue_func(n, p, k) < α # もしもpのP値が有意水準を下回ったら
                numtrials[i] = n # n回で目的を達成したことを記録に残す
                break
            end
        end
    end
    numtrials
end

# %%
numtrials_phack_exact = phacking(pvalue_exact, 10^3, 0.5)
[(n, ecdf(numtrials_phack_exact)(n)) for n in (10, 100, 1000)]

# %%
numtrials_phack_normal = phacking(pvalue_normal, 10^3, 0.5)
[(n, ecdf(numtrials_phack_normal)(n)) for n in (10, 100, 1000)]

# %%
"""ベイズ信頼区間ハックをプロット"""
function plot_bhack(; N=10^3, p=0.5, α = 0.05, L = 10^4, a = 1, b = a, kwargs...)
    numtrials_bhack = bayeshacking(N, p; α, L, a, b)
    numtrials_phack_exact = phacking(pvalue_exact, N, p; α, L)
    numtrials_phack_normal = phacking(pvalue_normal, N, p; α, L)

    P = plot(; legend=:bottomright)
    plot!(n -> ecdf(numtrials_bhack)(n), 0, N; label="Bayes hacking")
    plot!(n -> ecdf(numtrials_phack_exact)(n), 0, N; label="p-hacking (exact p-value)", ls=:dash)
    plot!(n -> ecdf(numtrials_phack_normal)(n), 0,N; label="p-hacking (normal dist. approx.)", ls=:dashdot)
    title!("ecdf(number of trials),  α = $α"; titlefontsize=12)
    plot!(; xlabel="n", ylabel="P(number of trials when stopped ≤ n)", guidefontsize=10)
    plot!(; kwargs...)

    Q = plot(Beta(a, b); ylim=(0, maximum(x -> pdf(Beta(a, b), x), 0.001:0.001:0.999) + 1), label="")
    plot!(; xtick=0:0.1:1, xlabel="parameter", ylabel="probability density")
    title!("prior = Beta($a, $b)"; titlefontsize=12)

    plot(P, Q; size=(640, 800), layout=(2, 1), leftmargin=2Plots.mm)
end

# %%
plot_bhack(; N=10^3, p=0.5, α = 0.05, a = 1, b = 1, xtick=0:100:1000)

# %%
plot_bhack(; N=10^3, p=0.5, α = 0.05, a = 0.5, b = 0.5, xtick=0:100:1000)

# %%
plot_bhack(; N=10^3, p=0.5, α = 0.05, a = 0.1, b = 0.1, xtick=0:100:1000)

# %%
plot_bhack(; N=10^3, p=0.5, α = 0.05, a = 10, b = 10, xtick=0:100:1000)

# %%
plot_bhack(; N=10^3, p=0.5, α = 0.05, a = 5, b = 15, xtick=0:100:1000)

# %% [markdown]
# 以下では通常のP値とベイズ版でのP値の類似物がほぼ等しくなることを確認しよう.

# %%
function pvalue_functions(; n = 20, k = 6, a = 1, b = a)
    m = k/n
    s = √(m*(1-m)/n)
    pmin, pmax = max(0.001, m - 4s), min(0.999, m + 4s)
    p = range(pmin, pmax; length=1000)

    pval_bayes = pvalue_bayes.(n, p, k; a, b)
    pval_exact = pvalue_exact.(n, p, k)
    pval_normal = pvalue_normal.(n, p, k)

    plot(;)
    plot!(p, pval_bayes; label="pvalue_bayes")
    plot!(p, pval_exact; label="pvalue_exact", ls=:dash)
    plot!(p, pval_normal; label="pvalue_normal", ls=:dashdot)
    title!("n = $n,  k = $k,  prior = Beta($a, $b)"; titlefontsize=12)
end

# %%
pvalue_functions(n = 10, k = 3)

# %%
pvalue_functions(n = 20, k = 6)

# %%
pvalue_functions(n = 50, k = 15)

# %%
pvalue_functions(n = 100, k = 30)

# %% [markdown]
# Stirlingの公式を使った計算によって, $k = np$ のとき, $n\to\infty$ で
#
# $$
# \begin{aligned}
# \text{log(marginal likelihood ratio)} &=
# \log B(k+a, n-k+b) - \log B(a, b)) - \log(p^k(1-p)^{n-p}
# \\ &=
# a\log p + b\log(1-p) - \frac{1}{2}\log\frac{np(1-p)}{2\pi} - \log B(a, b) + o(1)
# \end{aligned}
# $$
#
# 右辺は対数周辺尤度比の漸近挙動であり, $-(1/2)\log n + \mathrm{const.}$ の形をしている. 実際にはχ²分布のスケールに合わせるために２倍してあったことに注意せよ. そのおかげで, `k` が `Binomial(n, p)` 分布に従う確率変数のとき `logmarginallikrat(n, p, k; a, b)` はその最小値(上の２倍)と自由度１のχ²分布に従う確率変数の和に近似的に等しい.

# %%
"""
対数周辺尤度比 (χ²分布のスケールに合わせるために2倍しておく)
Bayes因子 = exp(logmarginallikrat(n, p, k; a, b)/2)
"""
@memoize function logmarginallikrat(n, p, k; a=1, b=a)
    logmarginallik = 2(logbeta(k + a, n - k + b) - logbeta(a, b))
    logmarginallik0 = 2(safemul(k, log(p)) + safemul(n - k, log(1 - p)))
    logmarginallik - logmarginallik0
end

"""`logmarginallikrat(n, p, k; a, b)` の最小値の定数項"""
lmlr_min_const(n, p; a=1, b=a) = 2(a*log(p) + b*log(1-p) - 1/2*log(p*(1-p)/(2π)) - logbeta(a, b))

"""`logmarginallikrat(n, p, k; a, b)` の最小値の対数項＋定数項"""
lmlr_min_approx(n, p; a=1, b=a) = -log(n) + lmlr_min_const(n, p; a, b)

"""`logmarginallikrat(n, p, k; a, b)` の分布のプロット"""
function plot_ecdflmlr(n, p=0.5; L=10^4, a=1, b=a)
    k = rand(Binomial(n, p), L)
    lmlr = @. logmarginallikrat(n, p, k; a, b)
    @show lmlr_min = logmarginallikrat.(n, p, n*p; a, b)
    @show lmlr_min_approx(n, p; a, b)
    @show -log(n)
    @show lmlr_min_const(n, p; a, b)
    lmlrmmin = lmlr .- lmlr_min
    plot(; legend=:bottomright)
    plot!(x -> ecdf(lmlrmmin)(x), 0, 8; label="ecdf(log marginal likelihood ratios - min)")
    plot!(x -> cdf(Chisq(1), x); label="cdf(Chisq(1), x)")
    title!("n = $n,  p = $p,  prior = Beta($a, $b)"; titlefontsize=12)
end

# %%
plot_ecdflmlr(10)

# %%
plot_ecdflmlr(100)

# %%
plot_ecdflmlr(1000)

# %%
"""
Bayes因子ハッキングに挑戦
Bayes因子 > threshold となるまでデータを取得し続ける。
"""
function try_bfhacking(N, p; threshold=10, L=10^4, a=1, b=a)
    numtrials = fill(N + 1, L) # N＋１で埋め尽くされた長さLのベクトル
    Threads.@threads for i in 1:L # 並列のでL回のシミュレーションを実行
        k = 0
        for n in 1:N # ベルヌイ試行回数nの最大値はN
            k += rand(Bernoulli(p)) # kに1または0を足す。1を足す確率はp.
            if logmarginallikrat(n, p, k; a, b) > 2log(threshold) # もしもベイズ因子が閾値を超えたら
                numtrials[i] = n # n回で目的を達成したこ"""logmarginallikrat(n, p, k; a, b)の最小値の定数項
                break
            end
        end
    end
    numtrials
end

"""Bayes因子ハッキングに類似のpハッキングに挑戦"""
function try_phacking_like_bfhacking(pvalue_func, N, p; threshold=10, L=10^4, a=1, b=a)
    numtrials = fill(N + 1, L) # N＋１で埋め尽くされた長さLのベクトル
    Threads.@threads for i in 1:L # 並列のでL回のシミュレーションを実行
        k = 0
        for n in 1:N # ベルヌイ試行回数nの最大値はN
            k += rand(Bernoulli(p)) # kに1または0を足す。1を足す確率はp.
            chisqthreshold_n = 2log(threshold) + log(n)
            α_n = ccdf(Chisq(1), chisqthreshold_n) # nごとに有意水準を設定する
            if pvalue_func(n, p, k) < α_n # もしもpのP値が有意水準を下回ったら
                numtrials[i] = n # n回で目的を達成したことを記録に残す
                break
            end
        end
    end
    numtrials
end

"""ベイズ因子ハッキングへの挑戦の様子をプロット"""
function plot_bfhacking(; N=5000, p=0.5, threshold=10, L=10^4, a=1, b=a)
    numtrials_bfhack = try_bfhacking(N, 0.5; threshold, L, a, b)
    numtrials_phack_like_bfhack_exact = try_phacking_like_bfhacking(pvalue_exact, N, 0.5; threshold, L, a, b)
    numtrials_phack_like_bfhack_normal = try_phacking_like_bfhacking(pvalue_normal, N, 0.5; threshold, L, a, b)

    P = plot(; legend=:bottomright)
    plot!(n -> ecdf(numtrials_bfhack)(n), 0, N; label="Bayes factor hacking")
    plot!(n -> ecdf(numtrials_phack_like_bfhack_exact)(n), 0, N; label="p-hacking like the above (exact)", ls=:dash)
    plot!(n -> ecdf(numtrials_phack_like_bfhack_normal)(n), 0, N; label="p-hacking like the above (normal dist. approx.)", ls=:dashdot)
    title!("ecdf(number of trials),  threshold = $threshold"; titlefontsize=12)
    plot!(; xlabel="n", ylabel="P(number of trials when stopped ≤ n)", guidefontsize=10)

    Q = plot(Beta(a, b); ylim=(0, maximum(x -> pdf(Beta(a, b), x), 0.001:0.001:0.999) + 1), label="")
    plot!(; xtick=0:0.1:1, xlabel="parameter", ylabel="probability density")
    title!("prior = Beta($a, $b)"; titlefontsize=12)

    plot(P, Q; size=(640, 800), layout=(2, 1), leftmargin=2Plots.mm)
end

# %% [markdown]
# threshold = 10 のとき、 Bayes因子 > 10 という条件でデータ取得を止めることに成功する確率は数パーセントのオーダーになっている。これの類似をP値を使ってやっても概ね結果は同じになる。有意水準を
#
# $$
# \alpha_n = \operatorname{ccdf}(\operatorname{Chisq}(1), 2\log(\mathrm{threshold}) + \log n)
# $$
#
# によって $n$ について単調減少するように決めている。

# %%
plot_bfhacking(; N=5000, p=0.5, threshold=10, a=1, b=1)

# %%
plot_bfhacking(; N=5000, p=0.5, threshold=10, a=0.5, b=0.5)

# %%
plot_bfhacking(; N=5000, p=0.5, threshold=10, a=0.1, b=0.1)

# %%
plot_bfhacking(; N=5000, p=0.5, threshold=10, a=10, b=10)

# %%
plot_bfhacking(; N=5000, p=0.5, threshold=10, a=5, b=15)

# %%

max_stars_repo_path: 0025/Bayes hacking.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0025/Bayes hacking.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0025/Bayes hacking.jl

max_forks_repo_name: genkuroki/public


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

JLLWrappers.@generate_main_file_header("polymake")
JLLWrappers.@generate_main_file("polymake", UUID("7c209550-9012-526c-9264-55ba7a78ba2c"))
end  # module polymake_jll


text:
export ft_ee_uawrite

function ft_ee_uawrite(ft_handle::Culong, buffer::Array{UInt8,1})
  len = length(buffer)
  ft_status = ccall((:FT_EE_UAWrite, d2xx),
                     Cuint,
                     (Culong, Ptr{UInt8}, Cuint),
                     ft_handle, buffer, len)
  checkstatus(ft_status)
  return nothing
end


text:
"""
    module Tessellations

a convenience wrapper to pacakge VoronoiDelaunay, with main export Tessellation.
"""
module Tessellations

export Tessellation,
tessellation_coords,
tessellation_inv_coords,
get_surrounding_sites

using ...Defaults
using RecipesBase
using VoronoiDelaunay

"""
    struct Tessellation

    Tessellation(x,y) -> t

create a Delaunay tessellation out of the points `x,y` (passed as vectors)
"""
struct Tessellation
    tessellation::DelaunayTessellation2D{Point2D}
    x::Array{Float64,1}
    y::Array{Float64,1}
    xy::Array{Point2D,1}
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64

    function Tessellation(x,y)
        xmin, xmax = extrema(x)
        xmin_new = (xmin+xmax)/2 - TESSELLATION_EXPANSION_FACTOR*(xmax-xmin)/2
        xmax_new = (xmin+xmax)/2 + TESSELLATION_EXPANSION_FACTOR*(xmax-xmin)/2
        ymin, ymax = extrema(y)
        ymin_new = (ymin+ymax)/2 - TESSELLATION_EXPANSION_FACTOR*(ymax-ymin)/2
        ymax_new = (ymin+ymax)/2 + TESSELLATION_EXPANSION_FACTOR*(ymax-ymin)/2
        return Tessellation(x,y,xmin_new,xmax_new,ymin_new,ymax_new)
    end
    function Tessellation(x,y,xmin,xmax,ymin,ymax)
        tessellation = DelaunayTessellation(length(x))
        xt, yt = tessellation_coords(x,y,xmin,xmax,ymin,ymax)
        p = Point.(xt,yt)
        push!(tessellation,deepcopy(p))
        return new(tessellation,x,y,p,xmin,xmax,ymin,ymax)
    end

    function Base.show(io::IO,t::Tessellation)
        Np = length(t.xy)
        Nt = length(t.tessellation._trigs)
        print(io,"Tessellation:")
        print(io," ",Np," points,")
        print(io," ",Nt," triangles")
    end
end

"""
    tessellation_coords(x,y,t::Tessellation) -> X,Y

map the point(s) `x,y` to the square ((1,2),(1,2)), needed for the VoronoiDelaunay pacakage.
"""
tessellation_coords(x,y,t::Tessellation) = tessellation_coords(x,y,t.xmin,t.xmax,t.ymin,t.ymax)
function tessellation_coords(x,y,args...)
    XT = Array{Float64}(undef,length(x))
    YT = Array{Float64}(undef,length(x))
    tessellation_coords!(XT,YT,x,y,args...)
    return XT, YT
end
function tessellation_coords!(XT,YT,x,y,args...)
    @inbounds for i ∈ eachindex(XT)
        XT[i],YT[i] = tessellation_coords(x[i],y[i],args...)
    end
    return nothing
end
function tessellation_coords(x::Real,y::Real,xmin,xmax,ymin,ymax)
    xt = min_coord + (max_coord-min_coord)*(x - xmin)/(xmax-xmin)
    yt = min_coord + (max_coord-min_coord)*(y - ymin)/(ymax-ymin)
    return xt, yt
end

"""
    tessellation_inv_coords(x,y,t::Tessellation) -> X,Y

undo the mapping to ((1,2),(1,2)) that is needed for the VoronoiDelaunay pacakage.
"""
tessellation_inv_coords(x,y,t::Tessellation) = tessellation_inv_coords(x,y,t.xmin,t.xmax,t.ymin,t.ymax)
function tessellation_inv_coords(xt,yt,xmin,xmax,ymin,ymax)
    X = Array{Float64}(undef,length(xt))
    Y = Array{Float64}(undef,length(xt))
    tessellation_inv_coords!(X,Y,xt,yt,xmin,xmax,ymin,ymax)
    return X, Y
end
function tessellation_inv_coords!(X,Y,xt,yt,xmin,xmax,ymin,ymax)
    @inbounds for i ∈ eachindex(X)
        X[i],Y[i] = tessellation_inv_coords(xt[i],yt[i],xmin,xmax,ymin,ymax)
    end
    return nothing
end
function tessellation_inv_coords(xt::Real,yt::Real,xmin,xmax,ymin,ymax)
    x = xmin + (xmax-xmin)*(xt - min_coord)/(max_coord-min_coord)
    y = ymin + (ymax-ymin)*(yt - min_coord)/(max_coord-min_coord)
    return x, y
end

"""
    get_surrounding_sites(t,x,y[,deep::Bool]) -> inds

`t` is a Tessellation, `x,y` is a positions.
returns indices of vertices in `t` that surround the point `x,y`

If `deep=true`, it goes 2-layers deep, meaning that it looks at the surrounding
triangles as well as those surrounding them.
"""
@inline function get_surrounding_sites(t::Tessellation,x::Real,y::Real,deep::Bool=false)
    t0 = locate(t,x,y)
    a, b, c = t0._a, t0._b, t0._c

    P = vcat(a,b,c)

    if isexternal(t0) || deep

        ta = movea(t,t0); Pa = [ta._a,ta._b,ta._c]
        tb = moveb(t,t0); Pb = [tb._a,tb._b,tb._c]
        tc = movec(t,t0); Pc = [tc._a,tc._b,tc._c]
        P = vcat(P,Pa,Pb,Pc)

        if deep
            taa = movea(t,ta); Paa = [taa._a,taa._b,taa._c]
            tab = moveb(t,ta); Pab = [tab._a,tab._b,tab._c]
            tac = movec(t,ta); Pac = [tac._a,tac._b,tac._c]

            tba = movea(t,tb); Pba = [tba._a,tba._b,tba._c]
            tbb = moveb(t,tb); Pbb = [tbb._a,tbb._b,tbb._c]
            tbc = movec(t,tb); Pbc = [tbc._a,tbc._b,tbc._c]

            tca = movea(t,tc); Pca = [tca._a,tca._b,tca._c]
            tcb = moveb(t,tc); Pcb = [tcb._a,tcb._b,tcb._c]
            tcc = movec(t,tc); Pcc = [tcc._a,tcc._b,tcc._c]

            P = vcat(P,Pa,Pb,Pc,Paa,Pab,Pac,Pba,Pbb,Pbc,Pca,Pcb,Pcc)
        end
        unique!(P)
    end

    inds = Array{Int}(undef,length(P))
    @inbounds for i ∈ eachindex(inds)
        temp = findfirst(isequal(P[i]),t.xy)
        inds[i] = isnothing(temp) ? 0 : temp
    end
    return inds[.!iszero.(inds)]
end

# extend VoronoiDelaunay methods to Tessellation
function VoronoiDelaunay.locate(t::Tessellation,x::Real,y::Real)
    @assert t.xmin ≤ x ≤ t.xmax && t.ymin ≤ y ≤ t.ymax "(x,y)=($x,$y) must be contained in ($(t.xmin),$(t.ymin)), ($(t.xmax),$(t.ymax))"
    return locate(t.tessellation,Point(tessellation_coords(x,y,t)...))
end
for fun ∈ (:movea, :moveb, :movec)
    @eval VoronoiDelaunay.$(fun)(tess::Tessellation,trig) = $(fun)(tess.tessellation,trig)
end

@recipe function f(t::Tessellation; delaunay=true, voronoi=false)
    aspect_ratio --> 1
    xlim --> (t.xmin,t.xmax)
    ylim --> (t.ymin,t.ymax)
    legend --> false

    if delaunay
        x, y = getplotxy(delaunayedges(t.tessellation))
        @series begin
            seriestype --> :path
            tessellation_inv_coords(x,y,t)
        end
    end
    if voronoi
        x, y = getplotxy(voronoiedges(t.tessellation))
        @series begin
            seriestype --> :path
            tessellation_inv_coords(x,y,t)
        end
    end
end

end


text:
VERSION >= v"0.6.0" && __precompile__(true)

module Iscalc

# Make sure you do
# using DataFrames
# using SQLite
# in that order in a fresh Julia session to ensure the conversion methods get defined

using Base.Dates, Requests, JuliaDB, CSV, FileIO, DataFrames, TextParse, HTTP, DataStreams, Missings, Roots, LaTeXStrings
using BusinessDays

#using IJulia, SQLite
#using Compat

include("GetData.jl")
include("OT.jl")
include("API.jl")

export obter, JCorridos, FluxosCaixa, fspot, fdesconto
export euronextBonds, factsheet, catalogo
export load, save, table, column, join, merge, stack, unstack, select, reduce   #JUliaDB
export latexstring # Latexstrings
println("\nISCALC V0.3.4 <2018-04-25>")

end # module

max_stars_repo_name: iscalprog/Iscalc

max_issues_repo_name: iscalprog/Iscalc

max_forks_repo_name: iscalprog/Iscalc


text:
#= -------------------------------------------------------------------------
# @file simulation_PID.jl
#
# @date 12/12/17 19:19:33
# @author Martin Noblia
# @email mnoblia@disroot.org
#
# @brief
#
# @detail
#
#  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#
module SimulationPID

#------------------------------------------------------------------------------
#                            parameters
#------------------------------------------------------------------------------
const  A      =  2.0
const  R      =  2.390057361376673e-5
const  cₑ     =  4.184e6
const  V      =  0.5
const  θ_ref  =  333.0
const  θ₀     =  281.0
const  θ_a    =  281.0
const  G_inv  = 0.4 * 836800.0
const  α      = cₑ * V
const  τ      = R * cₑ * V
const  q_nec = ((θ_ref - θ_a) / (R)) - A * G_inv

#------------------------------------------------------------------------------
#                            function system
#------------------------------------------------------------------------------

function update()
end

function update(system::Function, initial_condition::Number)
  @assert stop > start
  # Pre-allocate the timeseries matrix
  t = (start, stop)
  t_keep = collect(start:1.0:stop)

  # Pre-assign function
  θ_dot_comp(t, θ) = (G_inv * A / α) + (θ_a / τ) - (θ / τ) + (q_nec / α)

  # Perform the actual integration
  prob = ODEProblem(f, biomass, t)
  sol = solve(prob, saveat=t_keep, dense=false, save_timeseries=false)

end

function simulate(system::Function, start::Real, stop::Real, δ::Number)
   @assert stop > start
   output = allocate_simulation_out()
   while N > 1
      for time in start:δ:stop

      end
   end
end

function simulate(p, biomass; start::Int64=0, stop::Int64=500, use::Symbol=:stiff)
  @assert length(biomass) == size(p[:A],1)
  @assert use ∈ vec([:stiff :nonstiff])

  # Pre-allocate the timeseries matrix
  t = (float(start), float(stop))
  t_keep = collect(start:1.0:stop)

  # Pre-assign function
  f(t, y) = dBdt(t, y, p)

  # Perform the actual integration
  prob = ODEProblem(f, biomass, t)
  sol = solve(prob, saveat=t_keep, dense=false, save_timeseries=false, alg_hints=[use])

  output = Dict{Symbol,Any}(
  :p => p,
  :t => sol.t,
  :B => hcat(sol.u...)'
  )

  return output

end

end # end of SimulationPID


text:
using BenchmarkTools

const SUITE = BenchmarkGroup()
const PARAMS_FILE = joinpath(dirname(@__FILE__), "params.json")
const RESULTS_FILE = joinpath(dirname(@__FILE__), "results.json")

include(joinpath(dirname(@__FILE__), "MaxwellBenchmarks.jl"))
using .MaxwellBenchmarks
SUITE["Maxwell"] = MaxwellBenchmarks.SUITE

include(joinpath(dirname(@__FILE__), "VlasovBenchmarks.jl"))
using .VlasovBenchmarks
SUITE["Vlasov"] = VlasovBenchmarks.SUITE

if isfile(PARAMS_FILE)
    loadparams!(SUITE, BenchmarkTools.load(PARAMS_FILE)[1])
else
    tune!(SUITE)
    BenchmarkTools.save(PARAMS_FILE, params(SUITE))
end

r = run(SUITE["Vlasov"])
println(minimum(r))

if isfile(RESULTS_FILE)
    m = BenchmarkTools.load(RESULTS_FILE)[1]
    j = judge(minimum(r), minimum(m))
    println(j)
    if !isregression(j) && isimprovement(j)
        println("Saving new results")
        BenchmarkTools.save(RESULTS_FILE, r)
    end
else
    BenchmarkTools.save(RESULTS_FILE, r)
end


text:
function factored_norm(A::FacElem{NfOrdIdl, NfOrdIdlSet})
  b = Dict{fmpq, fmpz}()
  for (p, k) = A.fac
    n = norm(p)
    add_to_key!(b, fmpq(n), k)
    #if haskey(b, n)
    #  b[n] += k
    #else
    #  b[n] = k
    #end
  end
  bb = FacElem(b)
  simplify!(bb)
  return bb
end

function norm(A::FacElem{NfOrdIdl, NfOrdIdlSet})
  return evaluate(factored_norm(A))
end

function factored_norm(A::NfOrdFracIdl)
  n = norm(A)
  return FacElem(Dict(fmpq(numerator(n)) => 1, fmpq(denominator(n)) => -1))
end

function factored_norm(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  b = Dict{fmpq, fmpz}()
  for (p, k) = A.fac
    n = norm(p)
    v = numerator(n)
    add_to_key!(b, fmpq(v), k)
    #if haskey(b, v)
    #  b[v] += k
    #else
    #  b[v] = k
    #end
    v1 = denominator(n)
    add_to_key!(b, fmpq(v1), -k)
    #if haskey(b, v)
    #  b[v] -= k
    #else
    #  b[v] = -k
    #end
  end
  bb = FacElem(b)
  simplify!(bb)
  return bb
end

function norm(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  return evaluate(factored_norm(A))
end

const FacElemQ = Union{FacElem{fmpq, FlintRationalField}, FacElem{fmpz, FlintIntegerRing}}

function abs(A::FacElemQ)
  B = empty(A.fac)
  for (k,v) = A.fac
    ak = abs(k)
    add_to_key!(B, ak, v)
    #if haskey(B, ak)
    #  B[ak] += v
    #else
    #  B[ak] = v
    #end
  end
  if length(B) == 0
    return FacElem(Dict(one(base_ring(A)) => fmpz(1)))
  end
  return FacElem(B)
end

function ==(A::T, B::T) where {T <: FacElemQ}
  C = A*B^-1
  simplify!(C)
  return isone(C)
end

function isone(A::FacElemQ)
  C = simplify(A)
  return all(iszero, values(C.fac)) || all(isone, keys(C.fac))
end

function ==(A::NfOrdIdl, B::FacElem{NfOrdIdl, NfOrdIdlSet})
  C = inv(B)*A
  return isone(C)
end
==(B::FacElem{NfOrdIdl, NfOrdIdlSet}, A::NfOrdIdl) = A == B

function ==(A::NfOrdFracIdl, B::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  C = A*inv(B)
  return isone(C)
end
==(B::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, A::NfOrdIdl) = A == B

function isone(A::NfOrdFracIdl)
  B = simplify(A)
  return B.den == 1 && isone(B.num)
end

function ==(A::FacElem{NfOrdFracIdl,NfOrdFracIdlSet}, B::FacElem{NfOrdFracIdl,NfOrdFracIdlSet})
  return isone(A*inv(B))
end
function ==(A::FacElem{NfOrdIdl,NfOrdIdlSet}, B::FacElem{NfOrdIdl,NfOrdIdlSet})
  return isone(A*inv(B))
end
function ==(A::FacElem{NfOrdIdl,NfOrdIdlSet}, B::FacElem{NfOrdFracIdl,NfOrdFracIdlSet})
  return isone(A*inv(B))
end

==(A::FacElem{NfOrdFracIdl,NfOrdFracIdlSet}, B::FacElem{NfOrdIdl,NfOrdIdlSet}) = B==A

==(A::NfOrdFracIdl, B::FacElem{NfOrdIdl,NfOrdIdlSet}) = isone(A*inv(B))

function *(A::FacElem{NfOrdIdl,NfOrdIdlSet}, B::FacElem{NfOrdFracIdl,NfOrdFracIdlSet})
  C = deepcopy(B)
  for (i,k) = A.fac
    C *= FacElem(Dict(i//1 => k))
  end
  return C
end
*(A::FacElem{NfOrdFracIdl,NfOrdFracIdlSet}, B::FacElem{NfOrdIdl,NfOrdIdlSet}) = B*A

function isone(A::FacElem{NfOrdIdl, NfOrdIdlSet})
  A = simplify(A)
  return length(A.fac) == 1 && isone(first(keys(A.fac)))
end

function isone(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet})
  A = simplify(A)
  return length(A.fac) == 1 && (isone(first(keys(A.fac))) || iszero(first(values(A.fac))))
end


text:
export numberofstates

abstract type AbstractBalanceLaw{FT, D, S} end

Base.eltype(::AbstractBalanceLaw{FT}) where {FT} = FT
Base.ndims(::AbstractBalanceLaw{FT, D}) where {FT, D} = D
numberofstates(::AbstractBalanceLaw{FT, D, S}) where {FT, D, S} = S

function flux end
function wavespeed end
boundarystate(::AbstractBalanceLaw, n⃗, x⃗, q⁻, tag) = q⁻
source!(::AbstractBalanceLaw, dq, q, x⃗) = nothing
nonconservative_term!(::AbstractBalanceLaw, dq, q, x⃗) = nothing


text:
issimilar(a,b)=((typeof(a)==typeof(b)) && (size(a)==size(b)))

# Fix bug with deepcopy, where a shared bits array is copied multiple times:
# TODO: check if this bug is still there
Base.deepcopy_internal{T<:Number}(x::Array{T}, s::ObjectIdDict)=(haskey(s,x)||(s[x]=copy(x));s[x])

function Base.isapprox(x, y;
                       maxeps::Real = max(eps(eltype(x)), eps(eltype(y))),
                       rtol::Real=maxeps^(1/3), atol::Real=maxeps^(1/2))
    size(x) == size(y) || (warn("isapprox: $(size(x))!=$(size(y))"); return false)
    x = convert(Array, x)
    y = convert(Array, y)
    d = abs(x-y)
    s = abs(x)+abs(y)
    maximum(d - rtol * s) <= atol
end

# This is missing from Base
Base.convert{T,I}(::Type{Array{T,2}}, a::SparseMatrixCSC{T,I})=full(a)

# alternatives defined in cudart.jl
if !GPU
    typealias BaseArray{T,N} Union{Array{T,N},SubArray{T,N}}
end

copysync!(a::AbstractArray,b::AbstractArray)=copy!(a,b)
fillsync!(a::AbstractArray,x)=fill!(a,x)

# Define a more versatile version of randn!

import Base: randn!, GLOBAL_RNG

randn!{T}(A::AbstractArray{T}, mean=zero(T), std=one(T)) = randn!(GLOBAL_RNG, A, mean, std)

function randn!{T}(rng::AbstractRNG, A::AbstractArray{T}, mean=zero(T), std=one(T))
    for i in eachindex(A)
        @inbounds A[i] = mean+std*randn(rng)
    end
    A
end

# This is missing from sparse/linalg.jl: modified from (*) line 100.
import Base: A_mul_B!, A_mul_Bt!, At_mul_B!
using Base.LinAlg.BLAS: gemm!

### Add the ability to multiply arrays with other than 2 dimensions
mat2d(x)=(ndims(x)==2 ? x : (x2=reshape(x, size2(x));pointer(x2)===pointer(x)||error();x2))
A_mul_B!{T}(C::Array{T}, A::Array{T}, B::Array{T})=(gemm!('N','N',one(T),mat2d(A),mat2d(B),zero(T),mat2d(C)); C)
A_mul_Bt!{T}(C::Array{T}, A::Array{T}, B::Array{T})=(gemm!('N','T',one(T),mat2d(A),mat2d(B),zero(T),mat2d(C)); C)
At_mul_B!{T}(C::Array{T}, A::Array{T}, B::Array{T})=(gemm!('T','N',one(T),mat2d(A),mat2d(B),zero(T),mat2d(C)); C)

# y = w * x
function Base.A_mul_B!{TX,TvA,TiA}(Y::StridedMatrix{TX}, X::StridedMatrix{TX}, A::SparseMatrixCSC{TvA,TiA})
    mX, nX = size(X)
    nX == A.m || throw(DimensionMismatch())
    size(Y) == (mX, A.n) || throw(DimensionMismatch())
    fill!(Y,0)
    rowval = A.rowval
    nzval = A.nzval
    @inbounds for multivec_row=1:mX, col = 1:A.n, k=A.colptr[col]:(A.colptr[col+1]-1)
        Y[multivec_row, col] += X[multivec_row, rowval[k]] * nzval[k]
    end
    Y
end

# dw = dy * x'
function Base.A_mul_Bt!{TX,TvA,TiA}(Y::SparseMatrixCSC{TvA,TiA}, X::StridedMatrix{TX}, A::SparseMatrixCSC{TvA,TiA})
    error(:cpu_sparse_not_implemented_yet)
end

# dx = w' * dy is all dense

### DEAD CODE:

# # SIMILAR! create an array l.(n) similar to a given one.  If l.(n)
# # exists check and resize if necessary.

# function similar!(l, n, a, T=eltype(a), dims=size(a); fill=nothing)
#     if !isdefined(l,n) || (typeof(l.(n)) != typeof(a))
#         l.(n) = similar(a, T, dims)
#         fill != nothing && fillsync!(l.(n), fill)
#     elseif (size(l.(n)) != dims)
#         l.(n) = resize!(l.(n), dims)
#         fill != nothing && fillsync!(l.(n), fill)
#     end
#     return l.(n)
# end

# similar!(l, n, a, T, dims::Integer...; o...) = similar!(l,n,a,T,dims; o...)
# similar!(l, n, a, dims::Integer...; o...) = similar!(l,n,a,dims; o...)
# similar!(l, n, a, dims::Dims; o...) = similar!(l,n,a,eltype(a),dims; o...)

# issimilar1(a,b)=((eltype(a)==eltype(b)) && (isongpu(a)==isongpu(b)) && (length(a)==length(b)))
# issimilar2(a,b)=((eltype(a)==eltype(b)) && (isongpu(a)==isongpu(b)) && (size2(a)==size2(b)))

# This does not work in place!
# Base.resize!(a::Array, d::Dims)=similar(a, d)


text:
# ---
# title: 927. Three Equal Parts
# id: problem927
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Binary Search, Greedy
# link: <https://leetcode.com/problems/three-equal-parts/description/>
# hidden: true
# ---
#
# Given an array `A` of `0`s and `1`s, divide the array into 3 non-empty parts
# such that all of these parts represent the same binary value.
#
# If it is possible, return **any** `[i, j]` with `i+1 < j`, such that:
#
#   * `A[0], A[1], ..., A[i]` is the first part;
#   * `A[i+1], A[i+2], ..., A[j-1]` is the second part, and
#   * `A[j], A[j+1], ..., A[A.length - 1]` is the third part.
#   * All three parts have equal binary value.
#
# If it is not possible, return `[-1, -1]`.
#
# Note that the entire part is used when considering what binary value it
# represents.  For example, `[1,1,0]` represents `6` in decimal, not `3`.  Also,
# leading zeros are allowed, so `[0,1,1]` and `[1,1]` represent the same value.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,0,1,0,1]
#     Output: [0,3]
#
#
# **Example 2:**
#
#
#
#     Input: [1,1,0,1,1]
#     Output: [-1,-1]
#
#
#
# **Note:**
#
#   1. `3 <= A.length <= 30000`
#   2. `A[i] == 0` or `A[i] == 1`
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# This file was generated, do not modify it. # hide
using MLJ
using PyPlot
using PrettyPrinting
using StableRNGs
import DataFrames

MLJ.color_off() # hide
ioff() # hide
X, y = @load_boston
sch = schema(X)
p = length(sch.names)
n = sch.nrows
@show (n, p)
DataFrames.describe(y)


text:
module Optimise

export train!,
  SGD, ADAM, ADAMW, AdaMax, Momentum, Nesterov,
  RMSProp, ADAGrad, ADADelta, AMSGrad, NADAM, stop, StopException

struct Param{T}
  x::T
  Δ::T
end

Param(x::AbstractArray) = Param(x, zero(x))

include("optimisers.jl")
include("interface.jl")
include("train.jl")

using Flux.Tracker: TrackedArray

Param(x::TrackedArray) = Param(x.data, x.grad)
# Base.convert(::Type{Param}, x::TrackedArray) = Param(x.data, x.grad)

end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Tests related to ECEF to ECEF transformations using satellite state vector.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Get the current EOP Data.
#
# TODO: The EOP data obtained from IERS website does not match the values in the
# examples in [1]. However, this should be enough, because 1) the individual
# functions at the low level are tested using the same values of [1], and 2) the
# difference is smaller than 30 cm.

eop_iau1980  = read_iers_eop("./eop_IAU1980.txt",  Val(:IAU1980))
eop_iau2000a = read_iers_eop("./eop_IAU2000A.txt", Val(:IAU2000A))

# File: ./src/transformations/sv_ecef_to_ecef.jl
# ==============================================

# Functions: sv_ecef_to_ecef
# --------------------------

# ITRF <=> PEF
# ============

################################################################################
#                                  IAU-76/FK5
################################################################################

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
# one gets:
#
#   r_pef  = -1033.47503130   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_pef  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function sv_ecef_to_ecef ITRF <=> PEF" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    ## ITRF => PEF
    ## ===========

    r_itrf  = [-1033.4793830; 7901.2952754; 6380.3565958]
    v_itrf  = [-3.225636520; -2.872451450; +5.531924446]
    sv_itrf = orbsv(JD_UTC, r_itrf, v_itrf)
    sv_pef  = sv_ecef_to_ecef(sv_itrf, ITRF(), PEF(), JD_UTC, eop_iau1980)

    @test sv_pef.t === JD_UTC

    @test sv_pef.r[1] ≈ -1033.47503130 atol=1e-4
    @test sv_pef.r[2] ≈ +7901.30558560 atol=1e-4
    @test sv_pef.r[3] ≈ +6380.34453270 atol=1e-4

    @test sv_pef.v[1] ≈ -3.2256327470  atol=1e-7
    @test sv_pef.v[2] ≈ -2.8724425110  atol=1e-7
    @test sv_pef.v[3] ≈ +5.5319312880  atol=1e-7

    ## PEF => ITRF
    ## ===========

    r_pef   = [-1033.47503130; 7901.30558560; 6380.34453270]
    v_pef   = [-3.2256327470; -2.8724425110; +5.5319312880]
    sv_pef  = orbsv(JD_UTC, r_pef, v_pef)
    sv_itrf = sv_ecef_to_ecef(sv_pef, PEF(), ITRF(), JD_UTC, eop_iau1980)

    @test sv_itrf.t === JD_UTC

    @test sv_itrf.r[1] ≈ -1033.4793830 atol=1e-4
    @test sv_itrf.r[2] ≈ +7901.2952754 atol=1e-4
    @test sv_itrf.r[3] ≈ +6380.3565958 atol=1e-4

    @test sv_itrf.v[1] ≈ -3.225636520  atol=1e-7
    @test sv_itrf.v[2] ≈ -2.872451450  atol=1e-7
    @test sv_itrf.v[3] ≈ +5.531924446  atol=1e-7

end

################################################################################
#                                IAU-2006/2010
################################################################################

## ITRF <=> TIRS
## =============

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-14: Performing an IAU-2000 reduction [1, p. 220]
#
# According to this example and Table 3-6, using:
#
#   UTC    = April 6, 2004, 07:51:28.386009
#   r_itrf = -1033.4793830    i + 7901.2952754    j + 6380.3565958    k [km]
#   v_itrf =    -3.225636520  i -    2.872451450  j +    5.531924446  k [km/s]
#
# one gets the following:
#
#   r_tirs  = -1033.47503120   i + 7901.30558560   j + 6380.34453270   k [km]
#   v_tirs  =    -3.2256327470 i -    2.8724425110 j +    5.5319312880 k [km/s]
#
################################################################################

@testset "Function sv_ecef_to_ecef ITRF <=> TIRS" begin
    JD_UTC = date_to_jd(2004, 4, 6, 7, 51, 28.386009)

    ## ITRF => TIRS
    ## ===========

    r_itrf = [-1033.4793830; 7901.2952754; 6380.3565958]
    v_itrf = [-3.225636520; -2.872451450; +5.531924446]
    sv_itrf = orbsv(JD_UTC, r_itrf, v_itrf)
    sv_tirs = sv_ecef_to_ecef(sv_itrf, ITRF(), TIRS(), JD_UTC, eop_iau2000a)

    @test sv_tirs.t === JD_UTC

    @test sv_tirs.r[1] ≈ -1033.47503120 atol=1e-4
    @test sv_tirs.r[2] ≈ +7901.30558560 atol=1e-4
    @test sv_tirs.r[3] ≈ +6380.34453270 atol=1e-4

    @test sv_tirs.v[1] ≈ -3.2256327470  atol=1e-7
    @test sv_tirs.v[2] ≈ -2.8724425110  atol=1e-7
    @test sv_tirs.v[3] ≈ +5.5319312880  atol=1e-7

    ## TIRS => ITRF
    ## ===========

    r_tirs  = [-1033.47503120; 7901.30558560; 6380.34453270]
    v_tirs  = [-3.2256327470; -2.8724425110; +5.5319312880]
    sv_tirs = orbsv(JD_UTC, r_tirs, v_tirs)
    sv_itrf = sv_ecef_to_ecef(sv_tirs, TIRS(), ITRF(), JD_UTC, eop_iau2000a)

    @test sv_itrf.t === JD_UTC

    @test sv_itrf.r[1] ≈ -1033.4793830 atol=1e-4
    @test sv_itrf.r[2] ≈ +7901.2952754 atol=1e-4
    @test sv_itrf.r[3] ≈ +6380.3565958 atol=1e-4

    @test sv_itrf.v[1] ≈ -3.225636520  atol=1e-7
    @test sv_itrf.v[2] ≈ -2.872451450  atol=1e-7
    @test sv_itrf.v[3] ≈ +5.531924446  atol=1e-7

end


text:
using DAWG
using BenchmarkTools

#download csv dict from https://github.com/neologd/mecab-unidic-neologd/tree/master/seed
#decompress to any directory such as c:/tmp
#use only word with 3 to 11 characters (about 1.5 million words)

cd("C:/tmp")
sampledict = Dict()
f = open("mecab-unidic-user-dict-seed.20161215.csv")
for ln in eachline(f)
  word = (split(ln,","))[1]
  wlen = length(word)
  if wlen < 3 || wlen > 11
    continue
  end

  if !haskey(sampledict,word)
    sampledict[word]=true
  end
end
close(f)

l = collect(keys(sampledict))
println(length(l)," words in dictionary.")

@time smdawg = makeSparseMatrixDAWG(l)
#534.940754 seconds (673.87 M allocations: 30.877 GB, 18.57% gc time)

@time dadawg = makeDoubleArrayDAWG(l)

teststring = """
公開当日に生放送された。
出演者
木村拓哉
北川景子
陣内智則のプレゼンは時間の関係で省略になった（エンディングの浜辺のシーンで松たか子とのアドリブ）。
"""

smdawglist = DictionaryMatch(teststring,smdawg)
dadawglist = DictionaryMatch(teststring,dadawg)

@benchmark smdawglist = DictionaryMatch(teststring,smdawg,extract=false)

@benchmark dadawglist = DictionaryMatch(teststring,dadawg,extract=false)


text:
#ДАНО: На ограниченном внешней прямоугольной рамкой поле имеется ровно одна внутренняя перегородка в форме прямоугольника. Робот - в произвольной клетке поля между внешней и внутренней перегородками. РЕЗУЛЬТАТ: Робот - в исходном положении и по всему периметру внутренней перегородки поставлены маркеры.

function main6!(r::Robot)
    x, y1, y2 = getStartPos!(r)
    height = getHeight(r)
    up = moveToBorder!(r, x, height)
    mark_border!(r, up)
    setLeftDown!(r)
    while y2 > 0
          y2 -= 1
          move!(r, Nord)
    end
    while x > 0
          x -= 1
          move!(r, Ost)
    end
    while y1 > 0
          y1 -= 1
          move!(r, Nord)
    end
end

function getStartPos!(r::Robot)
    x, y1, y2 = 0, 0, 0
    while isborder(r, Sud) == false
          y1 += 1
          move!(r, Sud)
    end
    while isborder(r, West) == false
          x += 1
          move!(r, West)
    end
    while isborder(r, Sud) == false
          y2 += 1
          move!(r, Sud)
    end
    return x, y1, y2
end

function moveToBorder!(r::Robot, X, Y)
    x = 0
    while true
          y = 0
          if x % 2 == 0
                while isborder(r, Nord) == false
                      move!(r, Nord)
                      y += 1
                end
          else
                while isborder(r, Sud) == false
                      move!(r, Sud)
                      y += 1
                end
          end
          if y != Y || isborder(r, Ost) == true
                if x % 2 == 0
                      return true
                else
                      return false
                end
                break
          end
          move!(r, Ost)
          x += 1
    end
end

function getHeight(r::Robot)
    height = 0
    while isborder(r, Nord) == false
          height += 1
          move!(r, Nord)
    end
    while isborder(r, Sud) == false
          move!(r, Sud)
    end
    return height
end

function mark_border!(r, up)
    if up
          while isborder(r, Nord) == true
                putmarker!(r)
                move!(r, Ost)
          end
          move!(r, Nord)

          while isborder(r, West) == true
                putmarker!(r)
                move!(r, Nord)
          end
          move!(r, West)

          while isborder(r, Sud) == true
                putmarker!(r)
                move!(r, West)
          end
          move!(r, Sud)

          while isborder(r, Ost) == true
                putmarker!(r)
                move!(r, Sud)
          end
    else
          while isborder(r, Sud) == true
                putmarker!(r)
                move!(r, Ost)
          end
          move!(r, Sud)

          while isborder(r, West) == true
                putmarker!(r)
                move!(r, Sud)
          end
          move!(r, West)

          while isborder(r, Nord) == true
                putmarker!(r)
                move!(r, West)
          end
          move!(r, Nord)

          while isborder(r, Ost) == true
                putmarker!(r)
                move!(r, Nord)
          end
    end
end

function setLeftDown!(r::Robot)
    while isborder(r, Sud) == false
        move!(r, Sud)
    end
    while isborder(r, West) == false
        move!(r, West)
    end
end


text:
using Blink,RCall
w = Window()
loadurl(w, "file://"*pwd()*"/index.html")
sleep(1)

rfile = pwd()*"/gen.R"
@rput rfile
R"source(rfile)"

specs = ["CH4","CH3O2NO2","H2","N2O5","H2O2","CH3OH","HONO","CO","HO2NO2","CH3O","CH3NO3","HNO3","CH3OOH","HCHO","CH3O2","NO3","OH","NO2","NO","O3","HO2"]

data = []

for i in specs
    res = R"points($i,minno,maxno)"
    push!(data,[i,[res[1],res[2]]])

end

data

@js w data = $(data)
@js w JSON.stringify(data)

max_stars_repo_name: wolfiex/VisACC

max_issues_repo_name: wolfiex/VisACC

max_forks_repo_name: wolfiex/VisACC


text:
module OpenSpiel

using CxxWrap
using OpenSpiel_jll

import CxxWrap:argument_overloads

@wrapmodule(libspieljl)

include("patch.jl")

# export all
for n in names(@__MODULE__(); all=true)
    if Base.isidentifier(n) &&
        !startswith(String(n), "_") &&
        n ∉ (Symbol(@__MODULE__()), :eval, :include)
        @eval export $n
    end
end

function __init__()
    @initcxx
end

end # module


text:
using Atum
using Atum.EulerGravity

using CUDA
using LinearAlgebra: norm
using StaticArrays: SVector
using WriteVTK

import Atum: boundarystate
function boundarystate(law::EulerGravityLaw, n⃗, x⃗, q⁻, _)
  ρ⁻, ρu⃗⁻, ρe⁻ = EulerGravity.unpackstate(law, q⁻)
  ρ⁺, ρe⁺ = ρ⁻, ρe⁻
  ρu⃗⁺ = ρu⃗⁻ - 2 * (n⃗' * ρu⃗⁻) * n⃗
  SVector(ρ⁺, ρu⃗⁺..., ρe⁺)
end

function gravitywave(law, x⃗, t, add_perturbation=true)
  FT = eltype(law)
  x, z = x⃗

  Φ = grav(law) * z

  cv_d = FT(719)
  cp_d = γ(law) * cv_d
  R_d = cp_d - cv_d

  ΔT = FT(0.0001)
  H = FT(10e3)
  f = FT(0)
  L = FT(300e3)
  d = FT(5e3)
  x_c = FT(100e3)
  u_0 = FT(20)
  p_s = FT(1e5)
  T_ref = FT(250)

  g = grav(law)
  δ = g / (R_d * T_ref)
  c_s = sqrt(cp_d / cv_d * R_d * T_ref)
  ρ_s = p_s / (T_ref * R_d)

  if t == 0
    δT_b = ΔT * exp(-(x - x_c) ^ 2 / d ^ 2) * sin(π * z / H)
    δT = exp(δ * z / 2) * δT_b
    δρ_b = -ρ_s * δT_b / T_ref
    δρ = exp(-δ * z / 2) * δρ_b
    δu, δv, δw = 0, 0, 0
  else
    xp = x - u_0 * t

    δρ_b, δu_b, δv_b, δw_b, δp_b = zeros(SVector{5, Complex{FT}})
    for m in (-1, 1)
      for n in -100:100
        k_x = 2π * n / L
        k_z = π * m / H

        p_1 = c_s ^ 2 * (k_x ^ 2 + k_z ^ 2 + δ ^ 2 / 4) + f ^ 2
        q_1 = g * k_x ^ 2 * (c_s ^ 2 * δ - g) + c_s ^ 2 * f ^ 2 * (k_z ^ 2 + δ ^ 2 / 4)

        α = sqrt(p_1 / 2 - sqrt(p_1 ^ 2 / 4 - q_1))
        β = sqrt(p_1 / 2 + sqrt(p_1 ^ 2 / 4 - q_1))

        fac1 = 1 / (β ^ 2 - α ^ 2)
        L_m1 = (-cos(α * t) / α ^ 2 + cos(β * t) / β ^ 2) * fac1 + 1 / (α ^ 2 * β ^ 2)
        L_0 = (sin(α * t) / α - sin(β * t) / β) * fac1
        L_1 = (cos(α * t) - cos(β * t)) * fac1
        L_2 = (-α * sin(α * t) + β * sin(β * t)) * fac1
        L_3 = (-α ^ 2 * cos(α * t) + β ^ 2 * cos(β * t)) * fac1

        if α == 0
          L_m1 = (β ^ 2 * t ^ 2 - 1 + cos(β * t)) / β ^ 4
          L_0 = (β * t - sin(β * t)) / β ^ 3
        end

        δρ̃_b0 = -ρ_s / T_ref * ΔT / sqrt(π) * d / L *
                exp(-d ^ 2 * k_x ^ 2 / 4) * exp(-im * k_x * x_c) * k_z * H / 2im

        δρ̃_b = (L_3 + (p_1 + g * (im * k_z - δ / 2)) * L_1 +
              (c_s ^ 2 * (k_z ^ 2 + δ ^ 2 / 4) + g * (im * k_z - δ / 2)) * f ^ 2 * L_m1) * δρ̃_b0

        δp̃_b = -(g - c_s ^ 2 * (im * k_z + δ / 2)) * (L_1 + f ^ 2 * L_m1) * g * δρ̃_b0

        δũ_b = im * k_x * (g - c_s ^ 2 * (im * k_z + δ / 2)) * L_0 * g * δρ̃_b0 / ρ_s

        δṽ_b = -f * im * k_x * (g - c_s ^ 2 * (im * k_z + δ / 2)) * L_m1 * g * δρ̃_b0 / ρ_s

        δw̃_b = -(L_2 + (f ^ 2 + c_s ^ 2 * k_x ^ 2) * L_0) * g * δρ̃_b0 / ρ_s

        expfac = exp(im * (k_x * xp + k_z * z))

        δρ_b += δρ̃_b * expfac
        δp_b += δp̃_b * expfac

        δu_b += δũ_b * expfac
        δv_b += δṽ_b * expfac
        δw_b += δw̃_b * expfac
      end
    end

    δρ = exp(-δ * z / 2) * real(δρ_b)
    δp = exp(-δ * z / 2) * real(δp_b)

    δu = exp(δ * z / 2) * real(δu_b)
    δv = exp(δ * z / 2) * real(δv_b)
    δw = exp(δ * z / 2) * real(δw_b)

    δT_b = T_ref * (δp_b / p_s - δρ_b / ρ_s)
    δT = exp(δ * z / 2) * real(δT_b)
  end

  ρ_ref = ρ_s * exp(-δ * z)

  ρ = ρ_ref
  T = T_ref
  u = u_0
  w = FT(0)

  if add_perturbation
    ρ += δρ
    T += δT
    u += δu
    w += δw
  end

  e_kin = (u ^ 2 + w ^ 2) / 2
  e_int = cv_d * T
  ρe = ρ * (e_int + e_kin + Φ)

  return SVector(ρ, ρ * u, ρ * w, ρe)
end

function run(A, FT, N, KX, KY; outputvtk=true)
  Nq = N + 1

  law = EulerGravityLaw{FT, 2}()

  cell = LobattoCell{FT, A}(Nq, Nq)
  vx = range(FT(0), stop=FT(300e3), length=KX+1)
  vz = range(FT(0), stop=FT(10e3), length=KY+1)
  grid = brickgrid(cell, (vx, vz); periodic = (true, false))

  dg = DGSEM(; law, cell, grid, numericalflux = RoeFlux())
  #dg = ESDGSEM(; law, cell, grid,
  #             volume_numericalflux = EntropyConservativeFlux(),
  #             surface_numericalflux = RoeFlux())

  cfl = FT(1 // 3)
  dt = cfl * step(vz) / N / 330
  timeend = FT(30 * 60)

  q = gravitywave.(Ref(law), points(grid), FT(0))
  qref = gravitywave.(Ref(law), points(grid), FT(0), false)

  if outputvtk
    vtkdir = joinpath("output", "euler_gravity", "gravitywave")
    mkpath(vtkdir)
    pvd = paraview_collection(joinpath(vtkdir, "timesteps"))
  end

  do_output = function(step, time, q)
    if outputvtk && step % ceil(Int, timeend / 100 / dt) == 0
      filename = "step$(lpad(step, 6, '0'))"
      vtkfile = vtk_grid(joinpath(vtkdir, filename), grid)
      P = Bennu.toequallyspaced(cell)
      ρ, ρu, ρv, ρe = components(q)
      ρ_ref, ρu_ref, ρv_ref, ρe_ref = components(qref)
      vtkfile["δρ"] = vec(Array(P * (ρ - ρ_ref)))
      vtkfile["δρu"] = vec(Array(P * (ρu - ρu_ref)))
      vtkfile["δρv"] = vec(Array(P * (ρv - ρv_ref)))
      vtkfile["δρe"] = vec(Array(P * (ρe - ρe_ref)))
      vtk_save(vtkfile)
      pvd[time] = vtkfile
    end
  end

  odesolver = LSRK54(dg, q, dt)

  outputvtk && do_output(0, FT(0), q)
  solve!(q, timeend, odesolver; after_step=do_output)
  outputvtk && vtk_save(pvd)

  qexact = gravitywave.(Ref(law), points(grid), timeend)
  errf = map(components(q), components(qexact)) do f, fexact
    sqrt(sum(dg.MJ .* (f .- fexact) .^ 2))
  end
  norm(errf)
end

let
  A = CuArray
  FT = Float64
  N = 4

  ndof_x = 60
  ndof_y = 15
  KX_base = round(Int, ndof_x / N)
  KY_base = round(Int, ndof_y / N)

  nlevels = 2
  errors = zeros(FT, nlevels)
  for l in 1:nlevels
    KX = KX_base * 2 ^ (l - 1)
    KY = KY_base * 2 ^ (l - 1)
    errors[l] = run(A, FT, N, KX, KY)
  end
  rates = log2.(errors[1:(nlevels-1)] ./ errors[2:nlevels])
  @show errors
  @show rates
end


text:
# ---
# title: 1143. Longest Common Subsequence
# id: problem1143
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/longest-common-subsequence/description/>
# hidden: true
# ---
#
# Given two strings `text1` and `text2`, return the length of their longest
# common subsequence.
#
# A _subsequence_ of a string is a new string generated from the original string
# with some characters(can be none) deleted without changing the relative order
# of the remaining characters. (eg, "ace" is a subsequence of "abcde" while
# "aec" is not). A _common subsequence_  of two strings is a subsequence that is
# common to both strings.
#
#
#
# If there is no common subsequence, return 0.
#
#
#
# **Example 1:**
#
#
#
#     Input: text1 = "abcde", text2 = "ace"
#     Output: 3
#     Explanation: The longest common subsequence is "ace" and its length is 3.
#
#
# **Example 2:**
#
#
#
#     Input: text1 = "abc", text2 = "abc"
#     Output: 3
#     Explanation: The longest common subsequence is "abc" and its length is 3.
#
#
# **Example 3:**
#
#
#
#     Input: text1 = "abc", text2 = "def"
#     Output: 0
#     Explanation: There is no such common subsequence, so the result is 0.
#
#
#
#
# **Constraints:**
#
#   * `1 <= text1.length <= 1000`
#   * `1 <= text2.length <= 1000`
#   * The input strings consist of lowercase English characters only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@inline bitor(x,y) = Base.:(|)(x,y)
@inline bitand(x,y) = Base.:(&)(x,y)
@inline bitxor(x,y) = Base.:(⊻)(x,y)
@inline bitor(x,y,z) = bitor(x,bitor(y,z))
@inline bitand(x,y,z) = bitand(x,bitand(y,z))
@inline bitor(w,x,y,z) = bitor(bitor(w,x),bitor(y,z))
@inline bitand(w,x,y,z) = bitand(bitand(w,x),bitand(y,z))

# bitsof is like sizeof

bitsof(::Type{UInt8})   =   8
bitsof(::Type{UInt16})  =  16
bitsof(::Type{UInt32})  =  32
bitsof(::Type{UInt64})  =  64
bitsof(::Type{UInt128}) = 128
bitsof(::Type{Int8})    =   8
bitsof(::Type{Int16})   =  16
bitsof(::Type{Int32})   =  32
bitsof(::Type{Int64})   =  64
bitsof(::Type{Int128})  = 128

bitsof(x::U) where {U<:Unsigned} = bitsof(U)
bitsof(x::S) where {S<:Signed} = bitsof(S)

# with Markables, one bit is reserved for marking

bitsof(::Type{MarkUInt8})   =   7
bitsof(::Type{MarkUInt16})  =  15
bitsof(::Type{MarkUInt32})  =  31
bitsof(::Type{MarkUInt64})  =  63
bitsof(::Type{MarkUInt128}) = 127
bitsof(::Type{MarkInt8})    =   7
bitsof(::Type{MarkInt16})   =  15
bitsof(::Type{MarkInt32})   =  31
bitsof(::Type{MarkInt64})   =  63
bitsof(::Type{MarkInt128})  = 127

bitsof(x::U) where {U<:MarkableUnsigned} = bitsof(U)
bitsof(x::S) where {S<:MarkableSigned} = bitsof(S)

# lsbit(T) is used to mask the least significant bit of x::T

lsbit(::Type{UInt8})    = 0x01
lsbit(::Type{UInt16})   = 0x0001
lsbit(::Type{UInt32})   = 0x00000001
lsbit(::Type{UInt64})   = 0x0000000000000001
lsbit(::Type{UInt128})  = 0x00000000000000000000000000000001
lsbit(::Type{Int8})     = 0x01%Int8
lsbit(::Type{Int16})    = 0x0001%Int16
lsbit(::Type{Int32})    = 0x00000001%Int32
lsbit(::Type{Int64})    = 0x0000000000000001%Int64
lsbit(::Type{Int128})   = 0x00000000000000000000000000000001%Int128

lsbit(::Type{MarkUInt8})    = reinterpret(MarkUInt8, 0x01)
lsbit(::Type{MarkUInt16})   = reinterpret(MarkUInt16, 0x0001)
lsbit(::Type{MarkUInt32})   = reinterpret(MarkUInt32, 0x00000001)
lsbit(::Type{MarkUInt64})   = reinterpret(MarkUInt64, 0x0000000000000001)
lsbit(::Type{MarkUInt128})  = reinterpret(MarkUInt128, 0x00000000000000000000000000000001)
lsbit(::Type{MarkInt8})     = reinterpret(MarkInt8, 0x01)
lsbit(::Type{MarkInt16})    = reinterpret(MarkInt16, 0x0001)
lsbit(::Type{MarkInt32})    = reinterpret(MarkInt32, 0x00000001)
lsbit(::Type{MarkInt64})    = reinterpret(MarkInt64, 0x0000000000000001)
lsbit(::Type{MarkInt128})   = reinterpret(MarkInt128, 0x00000000000000000000000000000001)

# msbits(T) is used to mask all but the least significant bit of x::T

msbits(::Type{UInt8})    = 0xfe      # (-one(Int8))%UInt8 << 1
msbits(::Type{UInt16})   = 0xfffe
msbits(::Type{UInt32})   = 0xfffffffe
msbits(::Type{UInt64})   = 0xfffffffffffffffe
msbits(::Type{UInt128})  = 0xfffffffffffffffffffffffffffffffe
msbits(::Type{Int8})     = 0xfe%Int8
msbits(::Type{Int16})    = 0xfffe%Int16
msbits(::Type{Int32})    = 0xfffffffe%Int32
msbits(::Type{Int64})    = 0xfffffffffffffffe%Int64
msbits(::Type{Int128})   = 0xfffffffffffffffffffffffffffffffe%Int128

msbits(::Type{MarkUInt8})    = reinterpret(MarkUInt8, 0xfe)      # (-one(Int8))%UInt8 << 1
msbits(::Type{MarkUInt16})   = reinterpret(MarkUInt16, 0xfffe)
msbits(::Type{MarkUInt32})   = reinterpret(MarkUInt32, 0xfffffffe)
msbits(::Type{MarkUInt64})   = reinterpret(MarkUInt64, 0xfffffffffffffffe)
msbits(::Type{MarkUInt128})  = reinterpret(MarkUInt128, 0xfffffffffffffffffffffffffffffffe)
msbits(::Type{MarkInt8})     = reinterpret(MarkInt8, 0xfe)
msbits(::Type{MarkInt16})    = reinterpret(MarkInt16, 0xfffe)
msbits(::Type{MarkInt32})    = reinterpret(MarkInt32, 0xfffffffe)
msbits(::Type{MarkInt64})    = reinterpret(MarkInt64, 0xfffffffffffffffe)
msbits(::Type{MarkInt128})   = reinterpret(MarkInt128, 0xfffffffffffffffffffffffffffffffe)

# for isolation of mark, nomark bits
@inline lsbitof(x::T) where {T<:MarkableInteger}  =
    bitand(reinterpret(itype(T),x), lsbit(itype(T)))
@inline msbitsof(x::T) where {T<:MarkableInteger} =
    bitand(reinterpret(itype(T),x), msbits(itype(T)))
@inline lsbitof(x::T) where {T<:Union{Signed,Unsigned}}  =
    bitand(x, lsbit(itype(T)))
@inline msbitsof(x::T) where {T<:Union{Signed,Unsigned}} =
    bitand(x, msbits(itype(T)))

# for propogation of marks through operations
@inline any_lsbits(a::T) where {T<:MarkableInteger} = lsbitof(a)
@inline all_lsbits(a::T) where {T<:MarkableInteger} = lsbitof(a)
@inline any_lsbits(a::T, b::T) where {T<:MarkableInteger} = lsbitof(a) | lsbitof(b)
@inline all_lsbits(a::T, b::T) where {T<:MarkableInteger} = lsbitof(a) & lsbitof(b)
any_lsbits(a::T, b::T, c::T) where {T<:MarkableInteger} = bitor(lsbitof(a), lsbitof(b), lsbitof(c))
all_lsbits(a::T, b::T, c::T) where {T<:MarkableInteger} = bitand(lsbitof(a), lsbitof(b), lsbitof(c))
any_lsbits(a::T, b::T, c::T, d::T) where {T<:MarkableInteger} = bitand(lsbitof(a), lsbitof(b), lsbitof(c), lsbitof(d))
all_lsbits(a::T, b::T, c::T, d::T) where {T<:MarkableInteger} = bitand(lsbitof(a), lsbitof(b), lsbitof(c), lsbitof(d))

any_lsbits(a::T1, b::T2) where {T1<:MarkableInteger, T2<:MarkableInteger} = any_lsbits(promote(a,b)...,)
all_lsbits(a::T1, b::T2) where {T1<:MarkableInteger, T2<:MarkableInteger} = all_lsbits(promote(a,b)...,)


text:
function (+){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftSum2(a,b)
    eftRound(hi, lo, rounding)
end

function (-){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftDiff2(a,b)
    eftRound(hi, lo, rounding)
end

function (*){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftProd2(a,b)
    eftRound(hi, lo, rounding)
end

function (/){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftDiv2Approx(a,b)
    eftRound(hi, lo, rounding)
end

function (square){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSquare(a)
    eftRound(hi, lo, rounding)
end

function (sqrt){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSqrtApprox(a)
    eftRound(hi, lo, rounding)
end


text:
module SimulationParams

export SimulationParameters

struct SimulationParameters

    "Number of thermalization updates."
    burnin::Int

    "Number of simulation updates."
    nsteps::Int

    "Measurement frequncy."
    meas_freq::Int

    "Number of measurements made."
    num_meas::Int

    "Number of measurements averaged over in a bin"
    bin_size::Int

    "Total number of bins."
    num_bins::Int

    "Number of langevin steps per bin."
    bin_steps::Int

    "checkpointing frequency in seconds"
    chckpnt_freq::Int

    "path to where the data should be written"
    filepath::String

    "name of folder data will be dumped into"
    foldername::String

    "filepath + foldername"
    datafolder::String

    function SimulationParameters(burnin::Int, nsteps::Int, meas_freq::Int, num_bins::Int, chckpnt_freq::Real,
                                  filepath::String, foldername::String, datafolder::String)

        # sanity check
        @assert nsteps >= meas_freq * num_bins

        # calculating the number of measurements that will be made in the simulation
        @assert nsteps%max(1,meas_freq)==0
        @assert burnin%max(1,meas_freq)==0
        num_meas = div(nsteps, max(1,meas_freq) )

        # calculating the number of measurements that will be averaged over in each bin
        @assert num_meas%max(1,num_bins)==0
        bin_size = div(num_meas, max(1,num_bins) )

        # calculating the number of langevin time steps per bin
        bin_steps = meas_freq*bin_size

        # convert checkpointing frequency in minutes to seconds
        chckpnt_freq = round(Int,60*chckpnt_freq)

        new(burnin,nsteps,meas_freq,num_meas,bin_size,num_bins,bin_steps,chckpnt_freq,filepath,foldername,datafolder)
    end
end

end

max_stars_repo_name: cohensbw/Langevin

max_issues_repo_name: cohensbw/ElPhDynamics

max_forks_repo_name: cohensbw/ElPhDynamics