text:
#------------------------------------------------------------------------------

# returns f(x₀)
function evaluate_2(f::PolyElem, x₀)
    accum = convert(Int, lead(f))
    d = degree(f)
    reconstruct!(x₀)

    rec_x = rational_reconstruction(x₀)
    csr_x = sparse(to_dense(rec_x))

    I = sparse(to_dense(one(rec_x)))
    for i in 1 : d
        accum = accum*csr_x + convert(Int, coeff(f, d - i))*I
    end

    accum = from_dense(Array(accum), Nemo.QQ)

    return reconstruct!(modular_reduction(accum, base_ring(x₀)))
end

# returns f(x₀)
# O(d) multiplications/additions of x₀ if d = degree(f)
function evaluate(f::PolyElem, x₀)
    accum = lead(f)
    d = degree(f)
    reconstruct!(x₀)

    I = one(x₀)
    for i in 1 : d
        accum = accum * x₀ + coeff(f, d - i) * I
    end

    return reconstruct!(accum)
end

# returns f(x₀) * b
# expanding the brackets into:
#   f₀b + f₁x₀b + f₂x₀²b + ... + fₙx₀ⁿb
# and using the Horner scheme
function evaluate(f::PolyElem, x₀, b)
    E = one(x₀)
    accum = lead(f) * E * b
    d = degree(f)

    for i in 1 : d
        accum = x₀ * accum + coeff(f, d - i) * E * b
    end

    return accum
end

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

function shift_right_x(f::PolyElem)
    x = gen(parent(f))
    divexact(f - coeff(f, 0), x)
end

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

# Solves Ay = b
# Throws if A is singular
function square_nonsingular_randomized_wiedemann(A::AbstractSparseMatrix, b::AbstractSparseVector)
    # the function implements Algorithm 1 given on page 55
    # in https://doi.org/10.1109/TIT.1986.1057137
    # the notation and step numbering are taken from there

    field = base_ring(A)
    S, x = PolynomialRing(field, "x")
    n = order(A)

    # 1
    b₀ = b
    y = zero(b)
    k = 0
    d = 0

    subspace = [ b ]
    for i in 1 : 2n - 1
        push!(subspace, A * last(subspace))
    end

    # 2
    while !iszero(b)
        # 3
        u = random_sparsik(n, field, density=1)

        # 4
        seq = elem_type(field)[]
        for i in 1 : 2 * (n - d)
            push!(seq, inner(u, subspace[i]))
        end

        # 5
        f = minimal_polynomial(S(seq), x^(2 * (n - d)))

        if coeff(f, 0) == 0
            throw(SingularException(0))
        end

        f = f * inv(coeff(f, 0))

        # 6
        f⁻ = shift_right_x(f)

        accum = zero_sparsik(n, field)
        for j in 0 : degree(f⁻)
            reduce!(accum, subspace[j + 1], coeff(f⁻, j))
        end
        reduce!(y, accum, 1)
        b = b₀ + A * y
        d += degree(f)
        k += 1

        if d == n
            break
        end
    end

    return -y
end

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

# applies the polynomial `g` to the sequence `seq`
# by the definition from the 'Modern Computer Algebra'
#
# O(n^2) if n is the length of `seq`
# Actually, its O(nlogn): libflint implements FFT
function apply_polynomial(seq, g)
    S = parent(g)
    n = length(seq)
    d = degree(g)
    product = seq * reverse(g)
    S([coeff(product, i) for i in d : (2n - 1)])
end

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

function square_nonsingular_deterministic_wiedemann(A::AbstractSparseMatrix, b::AbstractSparseVector)
    # the function implements Algorithm 2 given on page 55
    # in https://doi.org/10.1109/TIT.1986.1057137
    # the notation and step numbering are taken from there

    field = base_ring(A)
    S, x = PolynomialRing(field, "x")
    n = order(A)

    # 1
    # O(nω) if ω is the number of nonzeroes in A
    subspace = [ b ]
    for i in 1 : 2n - 1
        push!(subspace, A * last(subspace))
    end

    # 2
    k = 0
    g = one(S)

    while k < n && degree(g) < n
        # 3
        u = unit_sparsik(n, k + 1, field)

        # 4
        # O(n)
        seq = [inner(x, u) for x in subspace]

        # 5, O(nlogn)
        d = degree(g)
        h = apply_polynomial(S(seq), g)

        # 6, O(n^2)
        f = minimal_polynomial(h, x^(2n - d))

        # 7, O(nlogn)
        g *= f

        # 8
        k += 1
    end

    if iszero(coeff(g, 0)) || isone(g)
        throw(SingularException(0))
    end

    # 9
    f = g * inv(coeff(g, 0))
    f⁻ = shift_right_x(f)

    accum = zero_sparsik(n, field)
    for j in 0 : degree(f⁻)
        reduce!(accum, subspace[j + 1], coeff(f⁻, j))
    end

    return -accum
end

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

# Returns a solution x to linear system Ax = b
# under an assumtion of A being square and nonsingular
# using Wiedemann system-solving algorithm as a base
# Throws if system is singular

"""
wiedemann_solve(A, b; proved)

Returns a solution x to linear system Ax = b under an assumtion of
A being square and nonsingular using Wiedemann system-solving algorithm.
By default, proved = true which calls 'square_nonsingular_deterministic_wiedemann'
otherwise calls 'square_nonsingular_randomized_wiedemann'.
"""
function wiedemann_solve(
        A::AbstractSparseMatrix,
        b::AbstractSparseVector;
        proved=true)

    if issquare(A)
        if proved
            return square_nonsingular_deterministic_wiedemann(A, b)
        else
            return square_nonsingular_randomized_wiedemann(A, b)
        end
    end

    error("nonsquare matrices are not supported")
end

# Returns the minimal annihilating polynomial of matrix A, that is,
# f ∈ F[x], such that
#   f(A)  = 0
#   f     = argmin(degree(f))
#
# The algorithm is randomized if subspace_minpoly is of W. family
# We should probably do some probabilistic analysis
# (or peek it at the W. paper)
function minimal_polynomial_wiedemann(
        A::AbstractSparseMatrix;
        subspace_minpoly=__deterministic_wiedemann_minpoly)

    S, _ = PolynomialRing(base_ring(A), "x")
    f = S(1)

    iterations = 4
    for _ in 1 : iterations
        f = lcm(f, subspace_minpoly(A, S))
    end

    return f
end

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

# Deterministic algorithm
# for minimal annihilating polynomial of complex matrix A
#
# Sadly, the source thyself does not contain any references
# about this Algorithm complexity analysis
#
#   O(nL(n))   ( O(n^4) basically )
# where n is the order of matrix
# and L(n) is the cost of multiplying two matrices of n×n
#
# Throws if the given matrix is zero
#
# `PolySpace` - a univariate polynomial ring to put the resulting polynomial to
function __deterministic_simple_minpoly(A, PolySpace)
    # the function implements the algorithm presented in paper
    #   An algorithm for the calculation of the minimal polynomial,
    #   S. BIAŁAS and M. BIAŁAS
    # the paper is available at
    #   http://bulletin.pan.pl/(56-4)391.pdf

    if iszero(A)
        throw(DomainError(A, "It is really zero"))
    end

    field = base_ring(A)
    n = order(A)

    # Computing the linear span of the powers of A
    columns = [ one(A) ]
    for i in 1 : n
        push!(columns, A * last(columns))
    end

    columns = Dict(i => vec(x) for (i, x) in enumerate(columns))

    B = from_rows(n+1, n^2, field, Array(1:n+1), columns)
    B = transpose!(B)

    V = linear_span!(collect(values(B.rows)))

    rows = Dict(i => x for (i, x) in enumerate(basis(V)))
    C = from_rows(length(rows), n + 1, field, Array(1:length(rows)), rows)
    C = transpose!(C)

    # Computing the degree of the minimal polynomial
    V = Subspacik(field)
    power = 0
    for i in 1:size(C, 1)
        eatcode = eat_sparsik!(V, deepcopy(C.rows[i]))
        if !haskey(C.rows, i) || eatcode == reduced
            power = i
            break
        end
    end

    # Computing the minimal polynomial
    x = gen(PolySpace)
    f = x^(power-1)
    last_row = haskey(C.rows, power) ? C.rows[power] : zero_sparsik(n + 1, field)
    for (piv, val) in last_row
        if haskey(V.echelon_form, piv)
            f -= val * x^(piv-1)
        end
    end

    return f
end

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

# Single deterministic iteration of minimal_polynomial algorithm
# Returns annihilating polynomial for {Aⁱb} for some random vector b
function __deterministic_wiedemann_minpoly(A::AbstractSparseMatrix, PolySpace)
    field = base_ring(A)
    n = order(A)
    b = random_sparsik(n, field, density=1.0)

    subspace = [ b ]
    for i in 1 : 2n - 1
        push!(subspace, A * last(subspace))
    end

    k = 0
    g = one(PolySpace)
    while k < n && degree(g) < n
        u = unit_sparsik(n, k + 1, field)

        # O(n)
        seq = [inner(x, u) for x in subspace]

        # O(nlogn)
        d = degree(g)
        h = apply_polynomial(PolySpace(seq), g)

        # O(n^2)
        f = minimal_polynomial(h, gen(PolySpace)^(2n - d))

        # O(nlogn)
        g *= f
        k += 1
    end

    return g
end

# Single randomized iteration of minimal_polynomial algorithm
# Returns annihilating polynomial for {Aⁱb} for some random vector b
function __randomized_wiedemann_minpoly(A::AbstractSparseMatrix, PolySpace)
    field = base_ring(A)
    n = order(A)
    b = random_sparsik(n, field, density=1.0)
    g = PolySpace(1)

    subspace = [ b ]
    for i in 1 : 2n - 1
        push!(subspace, A * last(subspace))
    end

    while degree(g) < n
        u = random_sparsik(n, field, density=1)

        seq = elem_type(field)[]
        for i in 1 : 2 * (n - degree(g))
            push!(seq, inner(u, subspace[i]))
        end

        f = PolySpace(0)
        try
            f = minimal_polynomial(
                PolySpace(seq),
                gen(PolySpace)^(2 * (n - degree(g)))
            )
        catch SingularException
            # singular sequence is expected here, so do not throw
            return g
        end

        g *= f
    end

    return g
end


text:
"""
    Lattice

A type defining a finite lattice in arbitary dimensions.

# Fields
$(TYPEDFIELDS)
"""
struct Lattice

    "Number of spatial dimensions."
    D::Int

    "Number of unit cells."
    N::Int

    "Linear extent of lattice in the direction of each lattice vector."
    L::Vector{Int}

    "Whether the lattice is periodic in the direction of each lattice vector."
    periodic::Vector{Bool}

    "Storage space for representing a location or displacement in the lattice."
    lvec::Vector{Int}
end

"""
    Lattice(L::Vector{Int}, periodic::Vector{Bool})

Constructs a [`Lattice`](@ref).
"""
function Lattice(L::Vector{Int},periodic::Vector{Bool})

    @assert 1 <= length(L) <= 3
    @assert length(L) == length(periodic)
    @assert all(l -> l > 0, L)

    # dimension of lattice
    D = length(L)
    # number of unit cells in lattice
    N = prod(L)
    # location/displacment array
    lvec = zeros(Int,D)

    return Lattice(D,N,L,periodic,lvec)
end

Lattice(; L, periodic) = Lattice(L,periodic)

"""
    Base.show(io::IO, lattice::Lattice)
    Base.show(io::IO, ::MIME"text/plain", lattice::Lattice)

Show lattice.
"""
Base.show(io::IO, lattice::Lattice) = print(io,"Lattice(D=$(lattice.D), N=$(lattice.N), ",
                                            lattice.L, ", ", lattice.periodic)
function Base.show(io::IO, ::MIME"text/plain", lattice::Lattice)

    (; D, N, L, periodic) = lattice
    println(io, "Lattice:")
    println(io, " • D = $D")
    println(io, " • N = $N")
    println(io, " • L = ", L)
    println(io, " • periodic = ", periodic)
    return nothing
end

"""
    valid_location(l::AbstractVector{Int}, lattice::Lattice)

Determine if `l` describes a valid location in the lattice.
"""
function valid_loc(l::AbstractVector{Int}, lattice::Lattice)

    (; D, N, L) = lattice

    isvalid = true
    # check if location valid in each direction
    for d in 1:D
        if !(0<=l[d]<L[d])
            isvalid = false
            break
        end
    end

    return isvalid
end

valid_loc(; l, lattice) = valid_loc(l, lattice)

"""
    pbc!(l::AbstractVector{Int}, lattice::Lattice)

Apply periodic boundary to unit cell location `l`.
"""
function pbc!(l::AbstractVector{Int}, lattice::Lattice)

    (; D, N, L, periodic) = lattice
    @assert length(l) == D
    for d in 1:D
        # check if given direction in lattice is periodic
        if periodic[d]
            # make sure each value is positive
            if l[d] < 0
                l[d] = l[d] + L[d] * (abs(l[d])÷L[d] + 1)
            end
            # apply periodic boundary conditions
            l[d] = mod(l[d], L[d])
        end
    end
    return nothing
end

pbc!(; l, lattice) = pbc!(l, lattice)

"""
    unitcell_to_loc!(l::AbstractVector{Int}, u::Int, lattice::Lattice)

Calculate the location `l` of a unit cell `u`.
"""
function unitcell_to_loc!(l::AbstractVector{Int}, u::Int, lattice::Lattice)

    (; D, N, L) = lattice
    @assert length(l) == D

    for d in D:-1:1
        N    = N ÷ L[d]
        l[d] = (u-1) ÷ N
        u    = mod1(u,N)
    end

    return nothing
end

unitcell_to_loc!(; l, u, lattice) = unitcell_to_loc!(l, u, lattice)

"""
    unitcell_to_loc(u::Int, lattice::Lattice)

Return the location of unit cell `u`.
"""
function unitcell_to_loc(u::Int, lattice::Lattice)

    l = zeros(Int,lattice.D)
    unitcell_to_loc!(l,u,lattice)
    return l
end

unitcell_to_loc(; u, lattice) = unitcell_to_loc(u, lattice)

"""
    loca_to_unitcell(l::AbstractVector{Int}, lattice::Lattice)

Return the unit cell found at location `l` in the lattice.
"""
function loc_to_unitcell(l::AbstractVector{Int}, lattice::Lattice)

    (; D, N, L) = lattice
    @assert length(l) == D
    @assert valid_loc(l, lattice)

    u = 1
    for d in D:-1:1
        N = N ÷ L[d]
        u = u + N * l[d]
    end

    return u
end

loc_to_unitcell(; l, lattice) = loc_to_unitcell(l, lattice)

ext: jl


text:
"""
    Scene TODO document this

## Constructors
$(SIGNATURES)

## Fields
$(FIELDS)
"""
mutable struct Scene <: AbstractScene
    "The parent of the Scene; if it is a top-level Scene, `parent == nothing`."
    parent

    "[`Events`](@ref) associated with the Scene."
    events::Events

    "The current pixel area of the Scene."
    px_area::Node{IRect2D}

    "Whether the scene should be cleared."
    clear::Bool

    "The `Camera` associated with the Scene."
    camera::Camera

    "The controls for the camera of the Scene."
    camera_controls::RefValue

    """
    The limits of the data plotted in this scene.
    Can't be set by user and is only used to store calculated data bounds.
    """
    data_limits::Node{Union{Nothing, FRect3D}}

    "The [`Transformation`](@ref) of the Scene."
    transformation::Transformation

    "The plots contained in the Scene."
    plots::Vector{AbstractPlot}
    # TODO why 2?
    theme::Attributes

    attributes::Attributes

    "Children of the Scene inherit its transformation."
    children::Vector{Scene}

    """
    The Screens which the Scene is displayed to.
    """
    current_screens::Vector{AbstractScreen}

    """
    Signal to indicate whether layouting should happen. If updated to true,
    the Scene will be layouted according to its attributes (`raw`, `center`, or `scale_plot`).
    """
    updated::Node{Bool}
end

function Scene(
        events::Events,
        px_area::Node{IRect2D},
        clear::Bool,
        camera::Camera,
        camera_controls::RefValue,
        scene_limits,
        transformation::Transformation,
        plots::Vector{AbstractPlot},
        theme::Attributes, # the default values a scene owns
        attributes::Attributes, # the actual attribute values of a scene
        children::Vector{Scene},
        current_screens::Vector{AbstractScreen},
        parent = nothing,
    )

    # indicates whether we can start updating the plot
    # will be set when displayed
    updated = Node(false)

    scene = Scene(
        parent, events, px_area, clear, camera, camera_controls,
        Node{Union{Nothing, FRect3D}}(scene_limits),
        transformation, plots, theme, attributes,
        children, current_screens, updated
    )

    onany(updated, px_area) do update, px_area
        if update && !(scene.camera_controls[] isa PixelCamera)
            if !scene.raw[]
                scene.update_limits[] && update_limits!(scene)
                scene.scale_plot[] && scale_scene!(scene)
                scene.center[] && center!(scene)
            end
        end
        nothing
    end
    if scene[:camera][] !== automatic && camera_controls[] == EmptyCamera()
        # camera shouldn't really be part of the attributes, especially since
        # it just adds the camera one time and after that isn't usable
        cam = pop!(scene.attributes, :camera)[]
        apply_camera!(scene, cam)
    end
    scene
end

function Scene(;clear = true, transform_func=identity, scene_attributes...)
    events = Events()
    theme = current_default_theme(; scene_attributes...)
    attributes = copy(theme)
    px_area = lift(attributes.resolution) do res
        IRect(0, 0, res)
    end
    on(events.window_area) do w_area
        if !any(x-> x ≈ 0.0, widths(w_area)) && px_area[] != w_area
            px_area[] = w_area
        end
    end
    scene = Scene(
        events,
        px_area,
        clear,
        Camera(px_area),
        RefValue{Any}(EmptyCamera()),
        nothing,
        Transformation(transform_func),
        AbstractPlot[],
        theme,
        attributes,
        Scene[],
        AbstractScreen[]
    )
    # Set the transformation parent
    scene.transformation.parent[] = scene
    current_global_scene[] = scene
    scene
end

function Scene(
        scene::Scene;
        events = scene.events,
        px_area = scene.px_area,
        clear = false,
        cam = scene.camera,
        camera_controls = scene.camera_controls,
        transformation = Transformation(scene),
        theme = copy(theme(scene)),
        current_screens = scene.current_screens,
        kw_args...
    )
    child = Scene(
        events,
        px_area,
        clear,
        cam,
        camera_controls,
        nothing,
        transformation,
        AbstractPlot[],
        merge(current_default_theme(), theme),
        merge!(Attributes(; kw_args...), scene.attributes),
        Scene[],
        current_screens,
        scene
    )
    push!(scene.children, child)
    child
end

function Scene(parent::Scene, area; clear=false, transform_func=identity, attributes...)
    events = parent.events
    px_area = lift(pixelarea(parent), convert(Node, area)) do p, a
        # make coordinates relative to parent
        IRect2D(minimum(p) .+ minimum(a), widths(a))
    end
    child = Scene(
        events,
        px_area,
        clear,
        Camera(px_area),
        RefValue{Any}(EmptyCamera()),
        nothing,
        Transformation(transform_func),
        AbstractPlot[],
        current_default_theme(; attributes...),
        merge!(Attributes(; attributes...), parent.attributes),
        Scene[],
        parent.current_screens,
        parent
    )
    push!(parent.children, child)
    return child
end

# Base overloads for Scene

Base.haskey(scene::Scene, key::Symbol) = haskey(scene.attributes, key)
Base.propertynames(scene::Scene) = fieldnames(Scene) ∪ propertynames(scene.attributes)

function Base.getindex(scene::Scene, key::Symbol)
    return haskey(scene.attributes, key) ? scene.attributes[key] : scene.theme[key]
end

function Base.setindex!(scene::Scene, value, key::Symbol)
    scene.attributes[key] = value
end

Base.parent(scene::Scene) = scene.parent
isroot(scene::Scene) = parent(scene) === nothing
function root(scene::Scene)
    while !isroot(scene)
        scene = parent(scene)
    end
    scene
end
parent_or_self(scene::Scene) = isroot(scene) ? scene : parent(scene)

Base.size(x::Scene) = pixelarea(x) |> to_value |> widths |> Tuple
Base.size(x::Scene, i) = size(x)[i]
function Base.resize!(scene::Scene, xy::Tuple{Number, Number})
    resize!(scene, IRect(0, 0, xy))
end
Base.resize!(scene::Scene, x::Number, y::Number) = resize!(scene, (x, y))
function Base.resize!(scene::Scene, rect::Rect2D)
    pixelarea(scene)[] = rect
end

"""
    getscreen(scene::Scene)
Gets the current screen a scene is associated with.
Returns nothing if not yet displayed on a screen.
"""
function getscreen(scene::Scene)
    if isempty(scene.current_screens)
        isroot(scene) && return nothing # stop search
        return getscreen(parent(scene)) # screen could be in parent
    end
    # TODO, when would we actually want to get a specific screen?
    return first(scene.current_screens)
end

getscreen(scene::SceneLike) = getscreen(rootparent(scene))

"""
    `update!(p::Scene)`

Updates a `Scene` and all its children.
Update will perform the following operations for every scene:
```julia
if !scene.raw[]
    scene.update_limits[] && update_limits!(scene)
    scene.scale_plot[] && scale_scene!(scene)
    scene.center[] && center!(scene)
end
```
"""
function update!(p::Scene)
    p.updated[] = true
    foreach(update!, p.children)
end

# Just indexing into a scene gets you plot 1, plot 2 etc
Base.iterate(scene::Scene, idx = 1) = idx <= length(scene) ? (scene[idx], idx + 1) : nothing
Base.length(scene::Scene) = length(scene.plots)
Base.lastindex(scene::Scene) = length(scene.plots)
getindex(scene::Scene, idx::Integer) = scene.plots[idx]
GeometryBasics.widths(scene::Scene) = widths(to_value(pixelarea(scene)))
struct Axis end

zero_origin(area) = IRect(0, 0, widths(area))

function child(scene::Scene; attributes...)
    Scene(scene, lift(zero_origin, pixelarea(scene)); attributes...)
end

"""
Creates a subscene with a pixel camera
"""
function cam2d(scene::Scene)
    return child(scene, clear = false, camera = cam2d!)
end

function campixel(scene::Scene)
    return child(scene, clear = false, camera = campixel!)
end

function getindex(scene::Scene, ::Type{Axis})
    for plot in scene
        isaxis(plot) && return plot
    end
    nothing
end

"""
Each argument can be named for a certain plot type `P`. Falls back to `arg1`, `arg2`, etc.
"""
function argument_names(plot::P) where P <: AbstractPlot
    argument_names(P, length(plot.converted))
end

function argument_names(::Type{<: AbstractPlot}, num_args::Integer)
    # this is called in the indexing function, so let's be a bit efficient
    ntuple(i-> Symbol("arg$i"), num_args)
end

function Base.empty!(scene::Scene)
    empty!(scene.plots)
    disconnect!(scene.camera)
    scene.data_limits[] = nothing
    scene.camera_controls[] = EmptyCamera()
    empty!(scene.theme)
    merge!(scene.theme, _current_default_theme)
    empty!(scene.children)
    empty!(scene.current_screens)
end

limits(scene::Scene) = scene.data_limits

limits(scene::SceneLike) = limits(parent(scene))

function scene_limits(scene::Scene)
    if scene.limits[] === automatic
        return scene.data_limits[]
    else
        return FRect3D(scene.limits[])
    end
end

# Since we can use Combined like a scene in some circumstances, we define this alias
theme(x::SceneLike, args...) = theme(x.parent, args...)
theme(x::Scene) = x.theme
theme(x::Scene, key) = x.theme[key]
theme(x::AbstractPlot, key) = x.attributes[key]
theme(::Nothing, key::Symbol) = current_default_theme()[key]

Base.push!(scene::Combined, subscene) = nothing # Combined plots add themselves uppon creation
function Base.push!(scene::Scene, plot::AbstractPlot)
    push!(scene.plots, plot)
    plot isa Combined || (plot.parent[] = scene)
    if !scene.raw[]
        # update scenes data limit for each new plot!
        scene.data_limits[] = if scene.data_limits[] === nothing
            data_limits(plot)
        else
            union(scene.data_limits[], data_limits(plot))
        end
    end
    for screen in scene.current_screens
        insert!(screen, scene, plot)
    end
end

function Base.delete!(screen::AbstractScreen, scene::Scene, plot::AbstractPlot)
    @warn "Deleting plots not implemented for backend: $(typeof(screen))"
end

function Base.delete!(scene::Scene, plot::AbstractPlot)
    len = length(scene.plots)
    filter!(x-> x !== plot, scene.plots)
    if length(scene.plots) == len
        error("$(typeof(plot)) not in scene!")
    end
    for screen in scene.current_screens
        delete!(screen, scene, plot)
    end
end

function Base.push!(scene::Scene, child::Scene)
    push!(scene.children, child)
    disconnect!(child.camera)
    nodes = map([:view, :projection, :projectionview, :resolution, :eyeposition]) do field
        lift(getfield(scene.camera, field)) do val
            push!(getfield(child.camera, field), val)
        end
    end
    cameracontrols!(child, nodes)
    return scene
end

events(scene::Scene) = scene.events
events(scene::SceneLike) = events(scene.parent)

camera(scene::Scene) = scene.camera
camera(scene::SceneLike) = camera(scene.parent)

cameracontrols(scene::Scene) = scene.camera_controls[]
cameracontrols(scene::SceneLike) = cameracontrols(scene.parent)

cameracontrols!(scene::Scene, cam) = (scene.camera_controls[] = cam)
cameracontrols!(scene::SceneLike, cam) = cameracontrols!(parent(scene), cam)

pixelarea(scene::Scene) = scene.px_area
pixelarea(scene::SceneLike) = pixelarea(scene.parent)

plots(scene::SceneLike) = scene.plots

const _forced_update_scheduled = Ref(false)

"""
Returns whether a scene needs to be updated
"""
function must_update()
    val = _forced_update_scheduled[]
    _forced_update_scheduled[] = false
    val
end

"""
Forces the scene to be re-rendered
"""
function force_update!()
    _forced_update_scheduled[] = true
end

const current_global_scene = Ref{Any}()

"""
Returns the current active scene (the last scene that got created)
"""
function current_scene()
    if isassigned(current_global_scene)
        current_global_scene[]
    else
        Scene()
    end
end

"""
Fetches all plots sharing the same camera
"""
plots_from_camera(scene::Scene) = plots_from_camera(scene, scene.camera)
function plots_from_camera(scene::Scene, camera::Camera, list = AbstractPlot[])
    append!(list, scene.plots)
    for child in scene.children
        child.camera == camera && !child.raw[] && plots_from_camera(child, camera, list)
    end
    list
end

"""
Flattens all the combined plots and returns a Vector of Atomic plots
"""
function flatten_combined(plots::Vector, flat = AbstractPlot[])
    for elem in plots
        if (elem isa Combined)
            flatten_combined(elem.plots, flat)
        else
            push!(flat, elem)
        end
    end
    flat
end

function insertplots!(screen::AbstractDisplay, scene::Scene)
    for elem in scene.plots
        insert!(screen, scene, elem)
    end
    foreach(s-> insertplots!(screen, s), scene.children)
end
update_cam!(scene::Scene, bb::AbstractCamera, rect) = nothing

function scale_scene!(scene::Scene)
    if is2d(scene) !== nothing && is2d(scene)
        area = pixelarea(scene)[]
        lims = scene_limits(scene)
        # not really sure how to scale 3D scenes in a reasonable way
        mini, maxi = minimum(lims), maximum(lims)
        l = ((mini[1], maxi[1]), (mini[2], maxi[2]))
        xyzfit = fit_ratio(area, l)
        s = to_ndim(Vec3f0, xyzfit, 1f0)
        scale!(scene, s)
    end
    return scene
end

function center!(scene::Scene, padding = 0.01)
    bb = boundingbox(scene)
    bb = transformationmatrix(scene)[] * bb
    w = widths(bb)
    padd = w .* padding
    bb = FRect3D(minimum(bb) .- padd, w .+ 2padd)
    update_cam!(scene, bb)
    scene
end
parent_scene(x::Combined) = parent_scene(parent(x))
parent_scene(x::Scene) = x

Base.isopen(x::SceneLike) = events(x).window_open[]

function is2d(scene::SceneLike)
    lims = scene_limits(scene)
    lims === nothing && return nothing
    return is2d(lims)
end
is2d(lims::Rect2D) = return true
is2d(lims::Rect3D) = widths(lims)[3] == 0.0

"""
    update_limits!(scene::Scene, limits::Union{Automatic, Rect} = scene.limits[], padding = scene.padding[])

This function updates the limits of the `Scene` passed to it based on its data.
If an actual limit is set by the theme or its attributes (scene.limits !== automatic),
it will not update the limits. Call update_limits!(scene, automatic) for that.
"""
update_limits!(scene::Scene) = update_limits!(scene, scene.limits[])

function update_limits!(scene::Scene, limits::Automatic, padding::Vec3f0 = scene.padding[])
    # for when scene is empty
    dlimits = data_limits(scene)
    dlimits === nothing && return #nothing to limit if there isn't anything
    tlims = (minimum(dlimits), maximum(dlimits))
    let avec = tlims[1], bvec = tlims[2]
        if !all(x-> all(isfinite, x), (avec, bvec))
            @warn "limits of scene contain non finite values: $(avec) .. $(bvec)"
            mini = map(x-> ifelse(isfinite(x), x, zero(x)), avec)
            maxi = Vec3f0(ntuple(3) do i
                x = bvec[i]
                ifelse(isfinite(x), x, avec[i] + oneunit(avec[i]))
            end)
            tlims = (mini, maxi)
        end
    end
    local new_widths
    let avec = tlims[1], bvec = tlims[2]
        new_widths = Vec3f0(ntuple(3) do i
            a = avec[i]; b = bvec[i]
            w = b - a
            # check for widths == 0.0... 3rd dimension is allowed to be 0 though.
            # TODO maybe we should allow any one dimension to be 0, and then use the other 2 as 2D
            with0 = (i != 3) && (w ≈ 0.0)
            with0 && @warn "Founds 0 width in scene limits: $(avec) .. $(bvec)"
            ifelse(with0, 1f0, w)
        end)
    end
    update_limits!(scene, FRect3D(tlims[1], new_widths), padding)
end

"""
    update_limits!(scene::Scene, new_limits::Rect, padding = Vec3f0(0))

This function updates the limits of the given `Scene` according to the given Rect.

A `Rect` is a generalization of a rectangle to n dimensions.  It contains two vectors.
The first vector defines the origin; the second defines the displacement of the vertices from the origin.
This second vector can be thought of in two dimensions as a vector of width (x-axis) and height (y-axis),
and in three dimensions as a vector of the width (x-axis), breadth (y-axis), and height (z-axis).

Such a `Rect` can be constructed using the `FRect` or `FRect3D` functions that are exported by
`AbstractPlotting.jl`.  See their documentation for more information.
"""
function update_limits!(scene::Scene, new_limits::Rect, padding::Vec3f0 = scene.padding[])
    lims = FRect3D(new_limits)
    lim_w = widths(lims)
    # use the smallest widths for scaling, to have a consistently wide padding for all sides
    minw = if lim_w[3] ≈ 0.0
        m = min(lim_w[1], lim_w[2])
        Vec3f0(m, m, 0.0)
    else
        Vec3f0(minimum(lim_w))
    end
    padd_abs = minw .* to_ndim(Vec3f0, padding, 0.0)
    scene.data_limits[] = FRect3D(minimum(lims) .- padd_abs, lim_w .+  2padd_abs)
    scene
end


text:
using Escher
using Gadfly

# See if using Gadfly initializes interop
@test length(methods(drawing)) > 2
@test method_exists(drawing, (Measures.Measure, Measures.Measure, Any))
@test method_exists(convert, (Type{Escher.Tile}, Gadfly.Plot))


text:
using SHA

@testset "Payload" begin

    host = "127.0.0.1"
    port = 8080
    server_url = "http://" * host * ":" * string(port)

    num_comm_rounds = 1
    fraction_clients = 0.1f0
    num_total_clients = 10

    base_config = Fed.Config.BaseConfig(
        server_url,
        num_comm_rounds,
        fraction_clients,
        num_total_clients
    )

    weights = rand(Float32, 199210)
    chunksize = 256
    is_client = false

    function serde(config)
        Fed.initialize_stats(config, length(weights))
        data = Fed.Serde.serialize_payload(config.payload_serde, weights)
        return Fed.Serde.deserialize_payload(config.payload_serde, data, "http://127.0.0.1:9090")
    end

    @testset "Vanilla" begin
        config = Fed.Config.VanillaConfig(base_config)
        @test serde(config) == weights
    end

    @testset "Quantized" begin
        config = Fed.Config.QuantizedConfig{UInt8}(base_config)
        res = serde(config)
        @test sum(res - weights) / length(weights) < 0.1 # quantzation generates a small error
    end

    @testset "Quantized diff" begin
        config = Fed.Config.QDiffConfig{UInt8}(base_config, chunksize, is_client)
        res = serde(config)
        @test sum(res - weights) / length(weights) < 0.1 # quantzation generates a small error
    end

    @testset "Quantized static diff" begin
        config = Fed.Config.QDiffStaticConfig{UInt8}(base_config, chunksize, is_client)
        res = serde(config)
        @test sum(res - weights) / length(weights) < 0.1 # quantzation generates a small error
    end

    @testset "GD" begin
        config = Fed.Config.GDConfig{UInt8}(base_config, chunksize, sha1, 0x05, host, 9090, is_client)
        res = serde(config)
        @test sum(res - weights) / length(weights) < 0.1 # quantzation generates a small error
    end
end


text:
"""
    shift_lattice(is)
Circular shift contiguous indices `is` in interval `[M÷2, M÷2)` such that set starts with 0,
where `M=length(is)`.

Inverse operation: [`shift_lattice_inv`](@ref).
"""
shift_lattice(is) = circshift(is, cld(length(is),2))

"""
    shift_lattice_inv(js)
Circular shift indices starting with 0 into a contiguous set in interval `[M÷2, M÷2)`,
where `M=length(js)`.

Inverse operation of [`shift_lattice`](@ref).
"""
shift_lattice_inv(js) = circshift(js, fld(length(js),2))

"""
    HubbardReal1DEP(address; u=1.0, t=1.0, v_ho=1.0)

Implements a one-dimensional Bose Hubbard chain in real space with external potential.

```math
\\hat{H} = -t \\sum_{\\langle i,j\\rangle} a_i^† a_j + \\sum_i ϵ_i n_i
+ \\frac{u}{2}\\sum_i n_i (n_i-1)
```

# Arguments

* `address`: the starting address, defines number of particles and sites.
* `u`: the interaction parameter.
* `t`: the hopping strength.
* `v_ho`: strength of the external harmonic oscillator potential ``ϵ_i = v_{ho} i^2``.
The first index is `i=0` and the maximum of the potential occurs in the centre of the
lattice.

# See also

* [`HubbardReal1D`](@ref)
* [`HubbardMom1D`](@ref)
* [`ExtendedHubbardReal1D`](@ref)

"""
struct HubbardReal1DEP{TT,A<:AbstractFockAddress,U,T,M} <: AbstractHamiltonian{TT}
    add::A
    ep::SVector{M,TT}
end

function HubbardReal1DEP(addr::BoseFS{<:Any,M}; u=1.0, t=1.0, v_ho=1.0) where M
    U, T, V = promote(float(u), float(t), float(v_ho))
    # js = range(1-cld(M,2); length=M)
    is = range(-fld(M,2); length=M) # [-M÷2, M÷2) including left boundary
    js = shift_lattice(is) # shifted such that js[1] = 0
    potential = SVector{M}(V*j^2 for j in js)
    return HubbardReal1DEP{typeof(U),typeof(addr),U,T,M}(addr, potential)
end

function Base.show(io::IO, h::HubbardReal1DEP)
    print(io, "HubbardReal1DEP($(h.add); u=$(h.u), t=$(h.t), v_ho=$(h.ep[2]))")
end

LOStructure(::Type{<:HubbardReal1DEP{<:Real}}) = IsHermitian()

Base.getproperty(h::HubbardReal1DEP, s::Symbol) = getproperty(h, Val(s))
Base.getproperty(h::HubbardReal1DEP{<:Any,<:Any,U}, ::Val{:u}) where U = U
Base.getproperty(h::HubbardReal1DEP{<:Any,<:Any,<:Any,T}, ::Val{:t}) where T = T
Base.getproperty(h::HubbardReal1DEP, ::Val{:add}) = getfield(h, :add)
Base.getproperty(h::HubbardReal1DEP, ::Val{:ep}) = getfield(h, :ep)

starting_address(h::HubbardReal1DEP) = h.add

function num_offdiagonals(::HubbardReal1DEP, address::BoseFS)
    return 2 * num_occupied_modes(address)
end

function diagonal_element(h::HubbardReal1DEP, address::BoseFS)
    sum(occupied_modes(address)) do index
        occnum, mode = index
        h.u * occnum * (occnum - 1) / 2 + h.ep[mode] * occnum
    end
end

function get_offdiagonal(h::HubbardReal1DEP, add::BoseFS, chosen)
    naddress, onproduct = hopnextneighbour(add, chosen)
    return naddress, - h.t * sqrt(onproduct)
end


text:
#FileIO2 code generation tools

#==Object builders
===============================================================================#

#==Function builders
===============================================================================#

#Generate a Base.show(File{}) Expr with a simplified (readable) signature.
function genexpr_ShowSimpleFileSign(fmt::Symbol)
	fmtstr = string(fmt)
	#TODO: how to escape fmt??? only know how to escape Expr
	return :(Base.show(io::IO, f::File{$fmt}) = print(io, "File{", $fmtstr, "}(\"$(f.path)\")"))
end

#==TODO:
Figure out how to make a macro in order to call genexpr_ShowSimpleFileSig using
nice syntax like:
	@Gen_ShowSimpleFileSign TextFmt HTMLFmt ...

#Generate Base.show(File{}) functions with a simplified (readable) signature.
macro Gen_ShowSimpleFileSig(fmtlist...)
	fnlist = Expr[]
	for fmt in fmtlist
		push!(fnlist, expr_FileShowSimple(fmt))
	end
	return [expr_FileShowSimple(fmt)) for fmt in fmtlist]
end
==#
#Last line


text:
using FinancialMonteCarlo, ForwardDiff;
@show "Black Scholes Model"
S0 = 100.0;
K = 100.0;
r = 0.02;
T = 1.0;
d = 0.01;
D = 90.0;

Nsim = 10000;
Nstep = 30;
sigma = 0.2
mc = MonteCarloConfiguration(Nsim, Nstep);
mc1 = MonteCarloConfiguration(Nsim, Nstep, FinancialMonteCarlo.AntitheticMC());
toll = 0.8

rfCurve = ZeroRate(r);

FwdData = Forward(T)
EUData = EuropeanOption(T, K)
AMData = AmericanOption(T, K)
BarrierData = BarrierOptionDownOut(T, K, D)
AsianFloatingStrikeData = AsianFloatingStrikeOption(T)
AsianFixedStrikeData = AsianFixedStrikeOption(T, K)
Model = BlackScholesProcess(sigma, Underlying(S0, d));

f(x::Vector) = pricer(BlackScholesProcess(x[1], Underlying(x[2], d)), ZeroRate(r), mc, EUData);

x = Float64[sigma, S0]

g = x -> ForwardDiff.gradient(f, x);

@show g(x)


text:
# -*- encoding: utf-8 -*-
#
# MIT License
#
# Copyright (c) 2016 Maurizio Tomasi
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
#
#
# # TERMINOLOGY
#
# The code in this file tries to follow the conventions laid in
# document LSPE-STRIP-SP-017 ("STRIP Telescope Coordinate Systems
# Specification") as much as possible. As the reader might not be
# acquainted with it, here is a short summary of the relevant
# terminology:
#
# +-----+-----------------------------+------------------------------------------+
# | SCS | Site Coordinate System      | Centered on the GPS station              |
# +-----+-----------------------------+------------------------------------------+
# | MCS | Mount Coordinate System     | Centered at the basis of the telescope   |
# +-------+---------------------------+------------------------------------------+
# | ACS | Azimuth Coordinate System   | Moveable, centered in the middle of      |
# |     |                             | ground motor. Positive angles go from    |
# |     |                             | North to *East* (beware!)                |
# +-----+-----------------------------+------------------------------------------+
# | ECS | Elevation Coordinate System | Moveable, the Y axis is aligned with     |
# |     |                             | the motor axis                           |
# +-----+-----------------------------+------------------------------------------+
# | BCS | Boresight Coordinate System | It should have been moveable, but the    |
# |     |                             | current design does not allow it (alas!) |
# +-----+-----------------------------+------------------------------------------+
# | RDP | Reference Detector Plane    | Centered in the focus of the telescope.  |
# |     |                             | All the beam vectors are specified in    |
# |     |                             | this frame.                              |
# +-----+-----------------------------+------------------------------------------+
#
# We always use a right-handed coordinate system:
#
#                 ^ Z axis (aligned with "up" direction, e.g., Zenith)
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#                 |
#          Origin O----------------> Y axis
#                /
#               /
#              /
#             /
#            /
#           v  X axis (aligned with North, when applicable)
#

import Healpix
import StaticArrays
import LinearAlgebra: ×, dot
import AstroLib
import Dates

export TENERIFE_LATITUDE_DEG, TENERIFE_LONGITUDE_DEG, TENERIFE_HEIGHT_M
export timetorotang, telescopetoground, groundtoearth
export genpointings!, genpointings, northdir, eastdir, polarizationangle

"Latitude of the LSPE/Strip site in Tenerife, in degrees"
const TENERIFE_LATITUDE_DEG = 28.30026

"Longitude of the LSPE/Strip site in Tenerife, in degrees"
const TENERIFE_LONGITUDE_DEG = -16.51012

"Height of the LSPE/Strip site in Tenerife, in meters"
const TENERIFE_HEIGHT_M = 2390

include("quaternions.jl")

"""
    timetorotang(time, rpm)

Convert a time into a rotation angle, given the number of rotations per minute.
The time should be expressed in seconds. The return value is in radians.
`time` can either be a scalar or a vector.
"""
function timetorotang(time_s, rpm)
    if rpm == 0
        0.0
    else
        2 * π * time_s * (rpm / 60)
    end
end

@doc raw"""
    telescopetoground(wheelanglesfn, time_s)

Return a quaternion of type `Quaternion{Float64}` representing the
coordinate transform from the focal plane to the ground of the
telescope. The parameter `wheelanglesfn` must be a function which
takes as input a time, `time_s`, in seconds, and it must return a
3-tuple containing the angles of the following motors, measured in
**radians**:

1. The boresight motor (rotation around the ``z`` axis,
   counterclockwise)

2. The altitude motor (rotation around the ``y`` axis,
   counterclockwise)

3. The ground motor (rotation around the ``z`` axis, **clockwise**:
   N→E→S→W)

# Example

`````julia
telescopetoground(3600.0) do
    # Boresight motor keeps a constant angle equal to 0°
    # Altitude motor remains at 20° from the Zenith
    # Ground motor spins at 1 RPM
    (0.0, deg2rad(20.0), timetorotang(time_s, 1))
end
`````
"""
function telescopetoground(wheelanglesfn, time_s)
    (wheel1ang, wheel2ang, wheel3ang) = wheelanglesfn(time_s)

    qwheel1 = qrotation_z(wheel1ang)
    qwheel2 = qrotation_y(wheel2ang)

    # The minus sign here takes into account the fact that the azimuth
    # motor requires positive angles to turn North into East
    qwheel3 = qrotation_z(-wheel3ang)

    qwheel3 * (qwheel2 * qwheel1)
end

"""
    groundtoearth(groundq, time_s, latitude_deg; day_duration_s=86400.0)

Return a quaternion of type `Quaternion{Float64}` representing the
coordinate transformation from the ground of the telescope to the
Equatorial coordinate system. The parameter `groundq` must be a
quaternion describing the coordinate transformation from the focal
plane of the telescope to the ground. The parameter `time_s` must be a
time in seconds, and `latitude_deg` is the latitude (in degrees, N is
positive) of the location where the observation is made.

The keyword `day_duration_s` specifies the length of a sidereal day in
seconds.
"""
function groundtoearth(groundq, time_s, latitude_deg; day_duration_s = 86400.0)
    locq = qrotation_y(deg2rad(90 - latitude_deg))
    earthq = qrotation_z(2π * time_s / day_duration_s)

    earthq * (locq * groundq)
end

"""
    vector2equatorial(dir, jd, latitude_deg, longitude_deg, height_m;
                      prec=true, nut=true, aber=true)

Transform the Healpix coordinates of a vector into equatorial coordinates.
The parameter `vector` is 3D vector, `jd` is the julian date. The
paramters `latitude_deg`, `longitude_deg` and `height_m` should contain the
latitude (in degrees, N is positive), the longitude (in degrees, counterclockwise
is positive) and the height (in meters) of the location where the observation is
made. The parameters `prec`, `nut` and `aber` allow to consider secondary effects
respectively of precession, nutation and aberration.
"""
function vector2equatorial(vector, jd, latitude_deg, longitude_deg, height_m,
                           prec, nut, aber, ref)
    (θ, ϕ) = Healpix.vec2ang(vector...)
    alt_rad = π / 2 - θ
    az_rad = 2π - ϕ

    ra_deg, dec_deg, _ = AstroLib.hor2eq(rad2deg(alt_rad),
                                         rad2deg(az_rad),
                                         jd,
                                         latitude_deg,
                                         longitude_deg,
                                         height_m,
                                         precession = prec,
                                         nutate = nut,
                                         aberration = aber,
                                         refract = ref)
    (deg2rad(dec_deg), deg2rad(ra_deg))
end

"""
    polarizationangle(northdir, eastdir, poldir)

Calculate the polarization angle projected in the sky in IAU
conventions. The parameters `northdir` and `eastdir` must be versors
that point towards the North and East, respectively; `poldir` must be
a versor that identify the polarization direction projected in the
sky. The return value is in radians, and it is zero if the
polarization angles points toward East, π/2 if it points toward North,
etc.

# Examples
```jldoctest
julia> polarizationangle([0, 0, 1], [0, 1, 0], [0, 1, 0])
0.0
julia> polarizationangle([0, 0, 1], [0, 1, 0], [0, 0, 1]) |> rad2deg
90.0
julia> polarizationangle([0, 0, 1], [0, 1, 0], [0, 0, -1]) |> rad2deg
-90.0
```
"""
function polarizationangle(northdir, eastdir, poldir)
    cosψ = clamp(dot(eastdir, poldir), -1, 1)
    sinψ = clamp(dot(northdir, poldir), -1, 1)
    atan(sinψ, cosψ)
end

northdir(θ, ϕ) = StaticArrays.@SArray [-cos(θ) * cos(ϕ), -cos(θ) * sin(ϕ), sin(θ)]
eastdir(θ, ϕ) = StaticArrays.@SArray [-sin(ϕ), cos(ϕ), 0.0]

"""
    northdir(θ, ϕ)
    eastdir(θ, ϕ)

Compute the North/East versor for a vector. The North for a vector v
is along the direction -dv/dθ, as θ is the colatitude and moves along
the meridian, and the East is along dv/dϕ.

# Examples
```jldoctest
julia> northdir(π/2, 0) ≈ [0, 0, 1]
true
julia> eastdir(π/2, 0) ≈ [0, 1, 0]
true
```
"""
northdir, eastdir

"""
    quat_to_angles(boreaxis, polaxis, quat)

Transform the boresight direction `dir` and the polarization direction
`poldir` according to quaternion `quat`. Return the 3-tuple (θ, ϕ, ψ)
representing the colatitude, longitude, and polarization angle
(calculated northward).

*Warning:* the definition of polarization angle does not work if the
observer is at the North or South Pole of the coordinate system. This
means that the polarization angle is not valid if the quaternion is
expressed in the ground reference system.

This function is used internally by [`genpointings`](@ref).
"""
function quat_to_angles(boreaxis, polaxis, quat)
    rotmatr = rotationmatrix_normalized(quat)
    boresight = rotmatr * boreaxis
    poldir = rotmatr * polaxis

    (θ, ϕ) = Healpix.vec2ang(boresight...)

    north = northdir(θ, ϕ)
    east = eastdir(θ, ϕ)
    (θ, ϕ, polarizationangle(north, east, poldir))
end

function genpointings!(wheelanglesfn,
                       beam_dir,
                       timerange_s,
                       dirs,
                       psi;
                       polaxis = Float64[1.0, 0.0, 0.0],
                       latitude_deg = TENERIFE_LATITUDE_DEG,
                       ground = false,
                       day_duration_s = 86400.0)

    if ground
        @assert size(dirs, 2) == 4
        @assert size(psi, 2) == 2
    else
        @assert size(dirs, 2) == 2
        @assert size(psi, 2) == 1
    end
    @assert size(dirs, 1) == size(psi, 1)

    for (idx, time_s) = enumerate(timerange_s)

        # This converts the RDP into the MCS (ground reference frame)
        groundq = telescopetoground(wheelanglesfn, time_s)
        # This converts the MCS into the celestial reference frame
        quat = groundtoearth(groundq, time_s, latitude_deg; day_duration_s = day_duration_s)
        θ, ϕ, curpsi = quat_to_angles(beam_dir, polaxis, quat)
        (dirs[idx, 1], dirs[idx, 2]) = (θ, ϕ)

        if ground
            # Re-run the transformation algorithm using the ground quaternion
            θ_ground, ϕ_ground, psi_ground = quat_to_angles(beam_dir, polaxis, groundq)

            (dirs[idx, 3], dirs[idx, 4]) = (θ_ground, ϕ_ground)
            (psi[idx, 1], psi[idx, 2]) = (curpsi, psi_ground)
        else
            psi[idx] = curpsi
        end
    end
end

function genpointings(wheelanglesfn,
                       beam_dir,
                       timerange_s;
                       polaxis = Float64[1.0, 0.0, 0.0],
                       latitude_deg = TENERIFE_LATITUDE_DEG,
                       ground = false,
                       day_duration_s = 86400.0)

    if ground
        dirs = Array{Float64}(undef, length(timerange_s), 4)
        psi = Array{Float64}(undef, length(timerange_s), 2)
    else
        dirs = Array{Float64}(undef, length(timerange_s), 2)
        psi = Array{Float64}(undef, length(timerange_s))
    end

    genpointings!(
        wheelanglesfn,
        beam_dir,
        timerange_s,
        dirs,
        psi;
        polaxis = polaxis,
        latitude_deg = latitude_deg,
        ground = ground,
        day_duration_s = day_duration_s,
    )

    (dirs, psi)
end

function genpointings!(wheelanglesfn,
                       beam_dir,
                       timerange_s,
                       t_start::Dates.DateTime,
                       skydirs,
                       skypsi;
                       polaxis = Float64[1.0, 0.0, 0.0],
                       latitude_deg = TENERIFE_LATITUDE_DEG,
                       longitude_deg = TENERIFE_LONGITUDE_DEG,
                       height_m = TENERIFE_HEIGHT_M,
                       precession = true,
                       nutation = true,
                       aberration = true,
                       refraction = true)

    @assert size(skydirs, 1) == size(skypsi, 1)
    @assert size(skydirs, 2) == 2
    @assert size(skypsi, 2) == 1

    for (idx, time_s) = enumerate(timerange_s)
        groundq = telescopetoground(wheelanglesfn, time_s)
        rotmatr = rotationmatrix_normalized(groundq)
        vector = rotmatr * beam_dir

        jd = AstroLib.jdcnv(t_start + Dates.Nanosecond(round(Int64, time_s * 1e9)))
        Dec_rad, Ra_rad = vector2equatorial(vector,
                                            jd,
                                            latitude_deg,
                                            longitude_deg,
                                            height_m,
                                            precession,
                                            nutation,
                                            aberration,
                                            refraction)

        poldir = rotmatr * polaxis
        north = northdir(π / 2 - Dec_rad, Ra_rad)
        east = eastdir(π / 2 - Dec_rad, Ra_rad)

        skydirs[idx, 1] = Dec_rad
        skydirs[idx, 2] = Ra_rad
        skypsi[idx] = polarizationangle(north, east, poldir)
    end
end

function genpointings!(wheelanglesfn,
                       beam_dir,
                       timerange_s,
                       t_start::Number,
                       skydirs,
                       skypsi;
                       kwargs...)
    genpointings!(
        wheelanglesfn,
        beam_dir,
        timerange_s,
        AstroLib.daycnv(t_start)::Dates.DateTime,
        skydirs,
        skypsi,
        kwargs...)
end

function genpointings(wheelanglesfn,
                      beam_dir,
                      timerange_s,
                      t_start::Dates.DateTime;
                      kwargs...)

    skydirs = Array{Float64}(undef, length(timerange_s), 2)
    skypsi = Array{Float64}(undef, length(timerange_s))

    genpointings!(
        wheelanglesfn,
        beam_dir,
        timerange_s,
        t_start,
        skydirs,
        skypsi;
        kwargs...)

    (skydirs, skypsi)
end

function genpointings(wheelanglesfn,
                       beam_dir,
                       timerange_s,
                       t_start::Number;
                       kwargs...)
    genpointings(
        wheelanglesfn,
        beam_dir,
        timerange_s,
        AstroLib.daycnv(t_start)::Dates.DateTime;
        kwargs...)
end

@doc raw"""
    genpointings!(wheelanglesfn, beam_dir, timerange_s, dirs, psi;
                  polaxis = Float64[1.0, 0.0, 0.0],
                  latitude_deg = TENERIFE_LATITUDE_DEG,
                  ground = false)
    genpointings(wheelanglesfn, beam_dir, timerange_s;
                 polaxis = Float64[1.0, 0.0, 0.0],
                 latitude_deg = TENERIFE_LATITUDE_DEG,
                 ground = false)
    genpointings!(wheelanglesfn, beam_dir, timerange_s, t_start, dirs, psi;
                  polaxis = Float64[1.0, 0.0, 0.0],
                  latitude_deg = TENERIFE_LATITUDE_DEG,
                  longitude_deg = TENERIFE_LONGITUDE_DEG,
                  height_m = TENERIFE_HEIGHT_M,
                  precession = true,
                  nutation = true,
                  aberration = true,
                  refraction = true)
    genpointings(wheelanglesfn, beam_dir, timerange_s, t_start;
                 polaxis=Float64[1.0, 0.0, 0.0],
                 latitude_deg=TENERIFE_LATITUDE_DEG,
                 longitude_deg=TENERIFE_LONGITUDE_DEG,
                 height_m=TENERIFE_HEIGHT_M,
                 precession = true,
                 nutation = true,
                 aberration = true,
                 refraction = true)

Generate a set of pointing directions for a STRIP detector. Each
function is provided in two flavours: the ones ending with `!` save
the results in the last two parameters `dirs` and `psi`, while the
others automatically allocate memory and return their results as a
pair `(dirs, psi)`.

The parameter `wheelanglesfn` must be a function which takes as input
a time in seconds and returns a 3-tuple containing the angles (in
radians) of the three motors:

1. The boresight motor (rotation around the ``z`` axis,
   counterclockwise)

2. The altitude motor (rotation around the ``y`` axis,
   counterclockwise)

3. The ground motor (rotation around the ``z`` axis, **clockwise**:
   N→E→S→W)

The meaning of the parameters/keywords is the following:

- `beam_dir` specifies the pointing direction of the mean (the boresight is
  [0, 0, 1]). It must be normalized.

- `timerange_s` is an enumerable type that specifies at which times
  (in seconds) pointings must be computed.

- `t_start` is a DateTime object or a Julian date, which specifies the
  UTC date and time when the observation starts

- `latitude_deg` is the latitude of the location where the observation
  is made (in degrees, North is positive). The default value is
  [`TENERIFE_LATITUDE_DEG`](@ref).

- `longitude_deg` is the longitude of the location where the
  observation is made. The default value is
  [`TENERIFE_LONGITUDE_DEG`](@ref).

- `height_m` is the elevation of the location where the observation is
  made (in meters). The default value is [`TENERIFE_HEIGHT_M`](@ref).

- If `ground` is `true`, the function will return a 4-tuple containing
  the colatitude and longitude measured in Equatorial coordinates
  (columns 1 and 2) and in ground coordinates (columns 3 and 4). If
  `ground` is `false`, only the Equatorial coordinates are computed.

- `polaxis` is the polarization axis; it must be normalized.

- `precession`: if `true` (the default), the Earth's precession is
  taken into account.

- `nutation`: if `true` (the default), the Earth's nutation is taken
  into account.

- `aberration`: if `true` (the default), stellar aberration is taken
  into account.

- `refraction`: if `true`, refraction corrections are taken into account.
  As these corrections are only valid for optical wavelengths, the
  default is `false`.

# Return values

If `t_start` is not provided, the function `genpointings` returns a
2-tuple containing the sky directions (a N×2 array containing
declination and right ascension, in Equatorial coordinates) and the
polarization angle for each time step. The function `genpointings!` sets
the values in the last two parameters `dirs` and `psi`.

If `t_start` is provided, the function `genpointings` returns a
2-tuple (4-tuple) containing the directions (a N×2 or Nx4 array
containing the declination and the right ascension) and the polarization
angles at each time step; `genpointings!` works as above.

# Examples

Here is an example using the form without `t_start`:

`````julia
dir, psi = genpointings([0, 0, 1], 0:0.1:1) do time_s
    # Boresight motor keeps a constant angle equal to 0°
    # Altitude motor remains at 20° from the Zenith
    # Ground motor spins at 1 RPM
    (0.0, deg2rad(20.0), timetorotang(time_s, 1))
end
`````

And here is an example using `t_start`; unlike the previous example,
we use a lambda function instead of the `do...end` notation.

`````julia
import Dates

dirs, psi = genpointings(time_s -> (0, deg2rad(20),
                                    timetorotang(time_s, 1)),
                         [0, 0, 1],
                         0:0.1:1,
                         Dates.DateTime(2022, 01, 01, 0, 0, 0),
                         latitude_deg=10.0,
                         longitude_deg=20.0,
                         height_m = 1000) do time_s
`````
"""
genpointings


text:
module TestTrajectory

include("preamble.jl")
include("toyenv.jl")
using LyceumBase: asvec, ntrajectories, nsamples, collate!

function makedata(e::AbstractEnvironment, lengths::Vector{Int})
    τs = [
        Trajectory(
            [rand(statespace(e)) for _ = 1:lengths[i]],
            [rand(observationspace(e)) for _ = 1:lengths[i]],
            [rand(actionspace(e)) for _ = 1:lengths[i]],
            rand(lengths[i]),
            rand(statespace(e)),
            rand(observationspace(e)),
            isodd(i),
        ) for i = 1:length(lengths)
    ]

    Bs = map(τs) do τ
        B = TrajectoryBuffer(e)
        nsamp = length(τ)
        resize!(B.S, nsamp)
        resize!(B.O, nsamp)
        resize!(B.A, nsamp)
        resize!(B.R, nsamp)
        resize!(B.sT, 1)
        resize!(B.oT, 1)
        resize!(B.done, 1)
        resize!(B.offsets, 2)

        B.S .= τ.S
        B.O .= τ.O
        B.A .= τ.A
        B.R .= τ.R
        B.sT[1] = τ.sT
        B.oT[1] = τ.oT
        B.done[1] = τ.done
        B.offsets[2] = length(τ)

        LyceumBase.checkrep(B)
    end

    return (Bs=Bs, τs=τs)
end

@testset "constructors" begin
    e = ToyEnv()
    @test_inferred TrajectoryBuffer(e)
    B = TrajectoryBuffer(e; sizehint = 123)
    @test length(B.S) == length(B.O) == length(B.A) == length(B.R) == 123
    @test length(B.sT) == length(B.oT) == length(B.done)
    @test length(B.offsets) - 1 == 0
end

@testset "collate" begin
    e = ToyEnv()
    let lens = [2,3], data = makedata(e, lens)
        B = TrajectoryBuffer(e)
        A = StructArray(collate!(B, data.Bs, sum(lens)))
        @test nsamples(B) == sum(lens)
        @test ntrajectories(B) == length(lens)

        @test length(flatten(A.S)) == length(B.S) == sum(lens)
        @test length(flatten(A.O)) == length(B.O) == sum(lens)
        @test length(flatten(A.A)) == length(B.A) == sum(lens)
        @test length(flatten(A.R)) == length(B.R) == sum(lens)
        @test length(A.sT) == length(lens)
        @test length(A.oT) == length(lens)
        @test length(A.done) == length(lens)

        @test all(1:length(A)) do i
            A[i].S == data.τs[i].S
        end
        @test all(1:length(A)) do i
            A[i].O == data.τs[i].O
        end
        @test all(1:length(A)) do i
            A[i].A == data.τs[i].A
        end
        @test all(1:length(A)) do i
            A[i].R == data.τs[i].R
        end
        @test all(1:length(A)) do i
            A[i].sT == data.τs[i].sT
        end
        @test all(1:length(A)) do i
            A[i].oT == data.τs[i].oT
        end
        @test all(1:length(A)) do i
            A[i].done == data.τs[i].done
        end
        @test B.offsets == [0, cumsum(lens)...]
    end
    let lens = [0,0], data = makedata(e, lens)
        B = TrajectoryBuffer(e, sizehint=123)
        collate!(B, data.Bs, 0)
        @test nsamples(B) == sum(lens)
        @test ntrajectories(B) == 0
        @test length(B.S) == 0
        @test length(B.O) == 0
        @test length(B.A) == 0
        @test length(B.R) == 0
        @test length(B.sT) == 0
        @test length(B.oT) == 0
        @test length(B.done) == 0
        @test B.offsets == [0]
    end
    let
        B = TrajectoryBuffer(e, sizehint=123)
        collate!(B, TrajectoryBuffer[], 0)
        @test nsamples(B) == 0
        @test ntrajectories(B) == 0
        @test length(B.S) == 0
        @test length(B.O) == 0
        @test length(B.A) == 0
        @test length(B.R) == 0
        @test length(B.sT) == 0
        @test length(B.oT) == 0
        @test length(B.done) == 0
        @test B.offsets == [0]
    end
    let lens = [0,0], data = makedata(e, lens)
        B = TrajectoryBuffer(e, sizehint=123)
        @test_throws ArgumentError collate!(B, data.Bs, -1)
        @test_throws ArgumentError collate!(B, data.Bs, 1)
    end
end

end # module


text:
######## REWRITE ################

push!(LOAD_PATH,"../src")

using Test
using TypesBasis
using ReadBasis

"""
		run_readbasis_test()

This is a function to test reading modified EMSL NWchem basis file.

Required File: ../examples/basissets/H.6-31G.mod

""" function run_readbasis_test()
	 @testset "Read Basis" begin
	     #H Gaussian 6-31
	     format="modnwchem";
	     filename="../examples/basissets/H.6-31G.mod";
	     basisref = [["     18.7311370              0.03349460       ",
	 	         "      2.8253937              0.23472695       ",
		         "      0.6401217              0.81375733       "],
		         ["      0.1612778              1.0000000        "]];
             nbasisfunc, basisfunc = readgaussparams(filename,format=format);

             @test size(basisref) == size(basisfunc)

             stringbasisref = filter(!isspace,join(vcat(basisref...)));
             stringbasisfunc = filter(!isspace,join(vcat(basisfunc...)));
             @test stringbasisref == stringbasisfunc;

	end #testset
end

"""
		run_readbasis_test()

This is a function to test the constructed gaussian basis datatype

Required File: ../examples/basissets/H.6-31G.mod

""" function run_retreadbasis_test()
    @testset "Return Read Basis" begin
        #H Gaussian 6-31
        basistype = "6-31"
        numprims = 4;
        format="modnwchem";
        filename="../examples/basissets/H.6-31G.mod"
        refexpo = [18.7311370,
                  2.8253937,
                  0.6401217,
                  0.1612778];
        refcoef = [0.03349460,
	 	   0.23472695,
		   0.81375733,
		   1.0000000];

        gaussbasis = retreadbasis(filename,format=format,
                                  btype=basistype)
        nbasisfunc = gaussbasis.nbasisfunc;
        basisfunc = gaussbasis.basisfunc;

        @test nbasisfunc == 2;

        @testset "Basis-Func. 1" begin
            @test basisfunc[1].info == basistype;
            @test basisfunc[1].nprims == 3;
            @test basisfunc[1].coefs == refcoef[1:3];
            @test basisfunc[1].alphas == refexpo[1:3];
         end

        @testset "Basis-Func. 2" begin
            @test basisfunc[2].info == basistype;
            @test basisfunc[2].nprims == 1;
            @test basisfunc[2].coefs[1] == refcoef[4];
            @test basisfunc[2].alphas[1] == refexpo[4];
        end

    end #testset
end #run_retreadbasis_test


text:
include("test_data.jl")

function test_cmp(io::Union{Nothing,IO} = nothing)
    map(TestData.test_cmp_data) do x
        ic = IntcodeMachine(x)
        run_intcode!(ic)
        res = fetch(ic)
        isnothing(res) && return @error "Intcode failed CMP '$x'"
        @info something(res)
        !isnothing(io) && flush(io)
    end
end

function test_jmp(io::Union{Nothing,IO} = nothing)
    map(TestData.test_jmp_data) do x
        ic = IntcodeMachine(x)
        run_intcode!(ic)
        res = fetch(ic)
        isnothing(res) && return @error "Intcode failed CMP '$x'"
        @info something(res)
        !isnothing(io) && flush(io)
    end
end

function test_larger(io::Union{Nothing,IO} = nothing)
    map(TestData.test_larger_data) do x
        ic = IntcodeMachine(x)
        run_intcode!(ic)
        res = fetch(ic)
        isnothing(res) && return @error "Intcode failed CMP '$x'"
        @info something(res)
        !isnothing(io) && flush(io)
    end
end


text:
using Requests, FTPClient, GZip, DataFrames

const WORKERS = 9 #seems like logical cores * (3/4) is the magic value
const CLEAN_WORKERS = true

if nworkers() < WORKERS || (nworkers()==1 && WORKERS == 1)
  addprocs(nworkers()==1?WORKERS-nworkers()+1:WORKERS-nworkers())
end

@everywhere using DataFrames

@everywhere if !isdefined(:constDef)
  const constDef = true
  const indexUrl = """http://www.sec.gov/Archives/edgar/full-index/"""
  const rootUrl = """http://www.sec.gov/Archives/"""
  const dataDir = "SEC10QK\\"
  const typeListToGet = ["10-Q","10-K"]
  const numMasterCols = 5
end

function getMaster(yearBegin::Int=2016, yearEnd::Int=2016)
  try mkdir(dataDir)  end

  #get the master file for each quarter
  for y::Int in yearBegin:yearEnd, q::Int in 1:4
    address::String = indexUrl * "$y" * "/QTR" * "$q" * "/master.gz"
    fileName = pwd() * "\\" * dataDir * "\\$y" * "QTR$q" * "master"
    println(fileName)
    try
      download(address, fileName * ".gz")
    catch err
      println(err)
    end

    #open and extract
    gHandle = GZip.open(fileName * ".gz")
    write(fileName, readlines(gHandle,chomp=false))

    close(gHandle)

  end

end

function parseMasterAsCSV2(yearBegin::Int=2016, yearEnd::Int=2016)

  #open and parse the master file for each quarter
  for y::Int in yearBegin:yearEnd, q::Int in 1:4

    fileName = pwd() * "\\" * dataDir * "\\$y" * "QTR$q" * "master"

    inFile::IOStream = open(fileName)
    outFile::IOStream = open(fileName*".txt", "w+")

    #get the number of lines to skip. #Hardcodes thes tructure
    for l in eachline(inFile)
      if ismatch(Regex(".txt"), l)
        lineArr::Vector{String} = split(l,'|')
        if (isnumber(lineArr[1]) == true && length(lineArr) == numMasterCols
          && sum(map(x::String->length(strip(x))>0?1:0, lineArr))==numMasterCols)
          write(outFile,join(lineArr,'\t'))
          write(outFile, '\n')
        end
      end
    end

    close(inFile)
    close(outFile)

  end
end

#helper function for pmap. Requires the year and quarter
@everywhere function _downloadFromMaster(y::Int, q::Int)
  mstrName::String = pwd() * "\\" * dataDir * "$y" * "QTR$q" * "master.txt"
  logFileName::String = "$(Dates.format(now(),"yymmdd_THMS"))_PID$(myid())_Log.txt"
  logFile::IOStream = open(logFileName, "w+")
  typeList::Vector{String} = deepcopy(typeListToGet) #make a local copy
  numFilesDLed = 0
  numErrors = 0
  numFiles = 0

  println(logFile, "Begin Processing ",mstrName)

  master::DataFrame = readtable(mstrName, header=false,separator='\t')
  println("\nSpun up process ", myid(),
    "\nDownloading from ", mstrName,
    "\nLog name: ", logFileName)
  for i in 1:(size(master,1)) #assumes first column is the CIK code
    if sum(map(x::String -> x==master[i,3]?1:0,typeList)) > 0
      outName::String = pwd() * "\\" * dataDir * "\\$y" * "QTR$q" *"_" * "$(master[i,1])" * ".txt"
      try
        download(rootUrl * "$(master[i,end])", outName)
        numFilesDLed += 1
      catch err
        println(logFile, "ERROR: ", err, "\nFailed to download ", outName)
        numErrors += 1
      end
      numFiles += 1
    end
  end
  println(logFile, "Successfully downloaded ", mstrName,
    "\n Success Count: ", numFilesDLed, "\nFailure Count:", numErrors,
    "\nTotal Files:", numFiles)
  println("Successfully downloaded ", mstrName,
    "\n Success Count: ", numFilesDLed, "\nFailure Count:", numErrors,
    "\nTotal Files:", numFiles)
end

function getSECFiles(yearBegin::Int=2016, yearEnd::Int=2016)

  #for each row and year
  yearByQuarter::Matrix{Int} = Matrix{Int}((yearEnd-yearBegin + 1)*4,2)
  for y::Int in yearBegin:yearEnd, q::Int in 1:4
    yearByQuarter[(y-yearBegin)*4+q, 1] = y
    yearByQuarter[(y-yearBegin)*4+q, 2] = q
  end

  pmap(_downloadFromMaster , yearByQuarter[:,1], yearByQuarter[:,2])
  println("All downloads attempted")

end

#getMaster()
#parseMasterAsCSV2()
getSECFiles()
cleanworkers()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
has_deps(x) = !isnothing(match(r"{(.*)}", x))

struct TypeDep
    data
    TypeDep(data) = new(replace(data, " " => ""))
end

function AbstractTrees.children(el::TypeDep)
    !has_deps(el.data) && return []

    dep = match(r"{(.*)}", el.data).captures[1]

    cb_count = 0
    siblings = []
    sibling_chars = Char[]
    for (i, char) ∈ enumerate(dep)
        cb_count += 1 * (char == '{')
        cb_count -= 1 * (char == '}')
        found = ((char == ',' || i == length(dep)) && cb_count == 0)
        if found
            if char == dep[end]
                push!(sibling_chars, char)
            end
            push!(siblings, String(sibling_chars))
            sibling_chars = Char[]
        else
            push!(sibling_chars, char)
        end
    end
    TypeDep.(replace.(siblings, r"<:\s*" => ""))
end

struct Converted
    initial_type
    final_type
end

extension_types = [
    "ANativeWindow",
    "AHardwareBuffer",
    "CAMetalLayer",
    "wl_surface",
    "wl_display",
    "Display",
    "VisualID",
    "xcb_connection_t",
    "xcb_window_t",
    "xcb_visualid_t",
    "MirConnection",
    "MirSurface",
    "HINSTANCE",
    "HWND",
    "RROutput",
]

type_conversions = Dict(
    "Cstring" => "Cstring",
    "Float32" => "Float32",
    "Float64" => "Float64",
    "Int32" => "Int32",
    "Int64" => "Int64",
    "Nothing" => "Nothing",
    "Ptr{Nothing}" => "Ptr{Nothing}",
    "UInt16" => "UInt16",
    "UInt32" => "UInt32",
    "UInt64" => "UInt64",
    "UInt8" => "UInt8",
    "Cdouble" => "Float64",
    "Cfloat" => "Float32",
    "Cint" => "Int32",
    "Csize_t" => "UInt",
    "Cssize_t" => "Int",
    "Cstring" => "String",
    "Cuint" => "UInt32",
    "VkBool32" => "Bool",
    "Ptr{Cstring}" => "Array{String}",
    "NTuple{256, UInt8}" => "String",
    "NTuple{16, UInt8}" => "String",
    "Ptr{Cvoid}" => "Ptr{Cvoid}",
    "Cvoid" => "Cvoid",
    "Ptr{ANativeWindow}" => "Ptr{Cvoid}",
    "Ptr{AHardwareBuffer}" => "Ptr{Cvoid}",
    "Ptr{CAMetalLayer}" => "Ptr{Cvoid}",
    "Ptr{xcb_connection_t}" => "Ptr{Cvoid}",
    "xcb_window_t" => "UInt32",
    "Window" => "UInt32",
    "Ptr{Display}" => "Ptr{Cvoid}",
    "void" => "Any"
    # "Ptr{Cfloat}" => "AbstractArray{<:Number}",
)

base_types = [
    "Bool",
    "String",
    "Function",
    "Int",
    "Int8",
    "Int16",
    "Int32",
    "Int64",
    "Cint",
    "UInt",
    "UInt8",
    "UInt16",
    "UInt32",
    "UInt64",
    "Cuint",
    "Csize_t",
    "Cssize_t",
    "Float32",
    "Float64",
    "Cfloat",
    "Cdouble",
    "Nothing",
    "Cstring",
    "Cvoid",
    "Number",
    "VersionNumber",
    "AbstractArray",
    "AbstractString",
    "Any",
    "",
]

function widen_type(type)
    stype = supertype(eval(Meta.parse(type)))
    stype <: Integer && return "Integer"
    stype <: AbstractFloat && return "Number"
    stype <: AbstractString && return "AbstractString"
    type
end

is_ptr_to_ptr(type) = startswith(type, "Ptr{Ptr{")
is_ptr(type) = startswith(type, "Ptr{")
is_ntuple(type) = startswith(type, "NTuple{")

is_vulkan_type(name) = any(startswith.(Ref(name), ["vk", "Vk", "VK_"]))
is_vulkan_struct(name) = name ∈ keys(api.structs)
is_extension_type(name) = name ∉ base_types && !is_vulkan_type(name)
is_extension_ptr(name) = is_ptr(name) && is_extension_type(inner_type(name))
is_base_type(name) = name ∈ base_types || is_literal(name)
is_literal(el) = occursin(r"[\.\"]+", el) || occursin(r"^\d+$", el)
is_expr(el) = occursin(r"[\(\)~\[\]]", el)
is_abstractarray_type(el) = startswith(el, "AbstractArray")
is_array_type(el) = startswith(el, "Array")
is_array_of_vk_objects(name, type, sname) = !isnothing(type) && is_ptr(type) && is_vulkan_struct(inner_type(type)) && is_array_variable(name, sname)

"Returns all the leaf types that a type depends on, including itself if the type contains no inner type."
function type_dependencies(type; include_type=false)
    deps = getproperty.(collect(Leaves(TypeDep(type))), :data)
end

"""
Take the inner type of a type, returning itself if the provided type has no dependencies.
"""
function inner_type(type)
    type_deps = type_dependencies(type, include_type=true)
    length(type_deps) == 1 && return last(type_deps)
    is_ntuple(type) && return last(type_deps) # ignore number of repetitions at index 1
    error("Cannot take inner type of $type")
end

is_opaque(ptr) = ptr == Ptr{Nothing}


text:
let
    K = initialize_regular_vtx_grid()
    G = initialize_grid_graph_from_vtx_grid(K)
    @test nv(G) == sum(K .!= 0)

    x_dim = size(K,1)
    y_dim = size(K,2)
    vtxs = [(i,j) for i in 1:x_dim for j in 1:y_dim if K[i,j] != 0]
    K2 = construct_vtx_grid(vtxs,(x_dim,y_dim))
    G2 = initialize_grid_graph_from_vtx_grid(K2)

    @test nv(G) == nv(G2)
    @test ne(G) == ne(G2)
end
let
    G = initialize_regular_grid_graph()
    initialize_grid_graph_with_obstacles((10,10),[(1,2),(3,4)])
    construct_vtx_grid([(1,2),(3,4)],(10,10))
end


text:
#

using Distributions
using Base.Test

dg = Gamma(1,1)

d_s = EdgeworthSum(dg,10)
d_m = EdgeworthMean(dg,10)
d_z = EdgeworthZ(dg,10)

dg_s = Gamma(10,1)
dg_m = Gamma(10,0.1)
dg_za = Gamma(10,1/std(dg_s))

for i = 0.01:0.01:0.99

    q = quantile(dg_s,i)
    @test isapprox(quantile(d_s, i)     , q          , atol=0.02)
    @test isapprox(cquantile(d_s, 1 - i), q          , atol=0.02)
    @test isapprox(cdf(d_s, q)          , i          , atol=0.002)
    @test isapprox(ccdf(d_s, q)         , 1 - i      , atol=0.002)
    @test isapprox(pdf(d_s, q)          , pdf(dg_s,q), atol=0.005)

    q = quantile(dg_m,i)
    @test isapprox(quantile(d_m, i)     , q          , atol=0.01)
    @test isapprox(cquantile(d_m, 1 - i), q          , atol=0.01)
    @test isapprox(cdf(d_m, q)          , i          , atol=0.002)
    @test isapprox(ccdf(d_m, q)         , 1 - i      , atol=0.002)
    @test isapprox(pdf(d_m, q)          , pdf(dg_m,q), atol=0.05)

    q = quantile(dg_za,i) - mean(dg_za)
    @test isapprox(quantile(d_z, i)     , q                       , atol=0.01)
    @test isapprox(cquantile(d_z, 1 - i), q                       , atol=0.01)
    @test isapprox(cdf(d_z, q)          , i                       , atol=0.002)
    @test isapprox(ccdf(d_z, q)         , 1 - i                   , atol=0.002)
    @test isapprox(pdf(d_z, q)          , pdf(dg_za,q+mean(dg_za)), atol=0.02)

end


text:
using AbInitioSoftwareBase: parentdir
using AbInitioSoftwareBase.Commands: MpiexecConfig
using AbInitioSoftwareBase.Inputs: Setter
using Dates: format, now
using EquationsOfStateOfSolids:
    EquationOfStateOfSolids, PressureEquation, Parameters, getparam, vsolve
using Express.EquationOfStateWorkflow: StOptim, ScfOrOptim
using QuantumESPRESSO.Commands: pw
using QuantumESPRESSO.Inputs.PWscf: PWInput, VerbositySetter, VolumeSetter, PressureSetter
using REPL.TerminalMenus: RadioMenu, request
using Setfield: @set!
using UnifiedPseudopotentialFormat  # To work with `download_potential`
using Unitful: Pressure, Volume, @u_str
using UnitfulAtomic

import Express.EquationOfStateWorkflow: MakeInput, FitEos, RunCmd

(::MakeInput{T})(template::PWInput, args...) where {T<:ScfOrOptim} =
    (customizer(args...) ∘ normalizer(T()))(template)

struct CalculationSetter <: Setter
    calc::ScfOrOptim
end
function (x::CalculationSetter)(template::PWInput)
    @set! template.control.calculation = if x.calc isa Scf  # Functions can be extended, not safe
        "scf"
    elseif x.calc isa StOptim
        "relax"
    else
        "vc-relax"
    end
    return template
end

normalizer(calc::ScfOrOptim) = VerbositySetter("high") ∘ CalculationSetter(calc)

struct OutdirSetter <: Setter
    timefmt::String
end
function (x::OutdirSetter)(template::PWInput)
    @set! template.control.outdir = abspath(
        joinpath(
            template.control.outdir,
            join((template.control.prefix, format(now(), x.timefmt), rand(UInt)), '_'),
        ),
    )
    mkpath(template.control.outdir)
    return template
end

customizer(volume::Volume, timefmt = "Y-m-d_H:M:S") =
    OutdirSetter(timefmt) ∘ VolumeSetter(volume)
function customizer(eos::PressureEquation, pressure::Pressure, timefmt = "Y-m-d_H:M:S")
    volumes = vsolve(eos, pressure)
    volume = length(volumes) > 1 ? _interactive_choose(volumes) : only(volumes)
    return OutdirSetter(timefmt) ∘ PressureSetter(pressure) ∘ VolumeSetter(volume)
end
customizer(params::Parameters, pressure::Pressure, timefmt = "Y-m-d_H:M:S") =
    customizer(PressureEquation(params), pressure, timefmt)

(x::RunCmd)(input, output = mktemp(parentdir(input))[1]; kwargs...) =
    pw(input, output; kwargs...)

function _interactive_choose(volumes)
    options = string.(volumes)
    menu = RadioMenu(options)
    choice = request("Choose the desired volume:", menu)
    choice == -1 ? throw(InterruptException()) : volumes[choice]
end


text:
const DEFAULT_FAVICON = "<link rel='icon' type='image/png' href='./favicon-16x16.png' sizes='16x16' /><link rel='icon' type='image/png' href='./favicon-32x32.png' sizes='32x32' />"
const DEFAULT_TITLE = "Swagger UI"

"""
    mutable struct Options

"""
mutable struct Options
    show_explorer::Bool
    custom_js::String
    custom_css::String
    custom_css_url::String
    custom_site_title::String
    custom_favicon::String
    swagger_options::Dict{String, Any}

    function Options(;
        show_explorer::Bool=true,
        custom_js::String="",
        custom_css::String="",
        custom_css_url::String="",
        custom_site_title::String=DEFAULT_TITLE,
        custom_favicon::String="",
        swagger_options::Dict{String, Any}=Dict{String, Any}())

        this = new()
        this.show_explorer = show_explorer
        this.custom_js = custom_js
        this.custom_css = custom_css
        this.custom_css_url = custom_css_url
        this.custom_site_title = custom_site_title
        this.custom_favicon = custom_favicon
        this.swagger_options = swagger_options
        return this
    end
end

function _options_to_json(options::Options, swagger_doc::Union{Dict{String, Any}, Nothing})::String
    swagger_options = options.swagger_options
    opts = Dict{String, Any}()
    opts["customOptions"] = swagger_options
    if !isnothing(swagger_doc)
        opts["swaggerDoc"] = swagger_doc
    end
    if haskey(swagger_options, "swaggerUrl")
        opts["swaggerUrl"] = swagger_options["swaggerUrl"]
    end
    if haskey(swagger_options, "swaggerUrls")
        opts["swaggerUrls"] = swagger_options["swaggerUrls"]
    end
    return JSON.json(opts)
end


text:
#---------------------------------------------------------------------
#*	2. conversion of a key component is known
#*--------------------------------------------------------------------
type stoic_conv_vap
	stoic_conv_vap()=begin
		new(
			stoic_vap(),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Key component",
				:Lower=>1,
				:Default=>1
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Molar conversion of key component",
				:Symbol=>"X_k"
			)),
			[
				:(_base_1.rate*_base_1._base_1._base_1.Tank.V = sumt(_base_1.stoic)/abs(sumt(_base_1.stoic(KComp,:)))*_base_1._base_1._base_1.Outletm.F*_base_1._base_1._base_1.Outletm.z(KComp)*kconv),
			],
			[
				"Reaction rate",
			],
			[:KComp,],
			[:kconv,]
		)
	end
	_base_1::stoic_vap
	KComp::DanaInteger
	kconv::DanaReal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export stoic_conv_vap
function setEquationFlow(in::stoic_conv_vap)
	addEquation(1)
end
function atributes(in::stoic_conv_vap,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/cstr"
	fields[:Brief]="Model of a generic vapour-phase stoichiometric CSTR based on conversion of a key component"
	fields[:Info]="
== Specify ==
* inlet stream
* conversion of a key component
"
	drive!(fields,_)
	return fields
end
stoic_conv_vap(_::Dict{Symbol,Any})=begin
	newModel=stoic_conv_vap()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
abstract TableRowDim
immutable HasRowDim <: TableRowDim end
immutable RowDimUnknown <: TableRowDim end

tblrowdim(tbl::AbstractTable) = RowDimUnknown()
nrow(tbl::AbstractTable) = _nrow(tbl, tblrowdim(tbl))
_nrow(tbl, ::RowDimUnknown) = error()
_nrow(tbl, ::HasRowDim) = ncol(tbl) > 0 ? length(columns(tbl)[1]) : 0

column_indexable(::AbstractTable) = false


text:
using Plots

f(z, C) = z * z + C

function convergence_judgment(C, initial_z)
    z = initial_z
    i::UInt8 = 0
    while i < 100
        z = f(z, C)
        if abs2(z) > 4
            return i
        end
        i += 1
    end
    return i
end

generate_grid(xrange, yrange, initial_z) =
    @. convergence_judgment(Complex(xrange', yrange), initial_z)

function main()
    N = 5000
    initial_z = 0.0
    xrange = yrange = range(-2, 2, length = N)
    @time grid = generate_grid(xrange, yrange, initial_z)
    heatmap(xrange, yrange, grid, aspect_ratio = 1, size = (500, 500))
end

main()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function initialize()
    #Initialize NO properties

    vvib_no = sqrt(e_vib_i*2.0/(mass_arr[1]*mass_arr[2]/(mass_arr[1] + mass_arr[2])))		# Initial vibrational velocity of NO
    vz_no = sqrt(e_trans_i*2.0/(mass_arr[1] + mass_arr[1]))				# Initial z velocity of NO

    #initial position of NO
    x0no = zeros(Float64, 2, 3)

    xno = rand() #temporary variable
    xno = @. xno * [cell[1], cell[2], 1.0]
    xno[3] = xno[3] + 10.5 # Initial z position of NO = 10

    theta_no = rand() # Initial orientation of NO, 0 = O-down, PI = N-down
    theta_no = asin(2.0 * theta_no - 1.0) + π/2.0
    phi_no = rand()
    phi_no = 2.0*π*phi_no

    x0no[1, :] = xno + r_no/2.0 * [sin(phi_no)*sin(theta_no), cos(phi_no)*sin(theta_no), cos(theta_no)]
    x0no[2, :] = xno - r_no/2.0 * [sin(phi_no)*sin(theta_no), cos(phi_no)*sin(theta_no), cos(theta_no)]

    xno = x0no[1,:] - x0no[2, :]
    xno = xno / sqrt(xno'*xno)

    #initial velocity of no
    v0no = zeros(Float64, 2, 3)
    phi_no = 0.0
    v0no[1, :] = vz_no * [sin(phi_no)*sin(phi_inc), -cos(phi_no)*sin(phi_inc), -cos(phi_inc)]
    v0no[2, :] = v0no[1, :]
    v0no[1, :] = v0no[1, :] + vvib_no * xno * mass_arr[2] / (mass_arr[1] + mass_arr[2])
    v0no[2, :] = v0no[2, :] - vvib_no * xno * mass_arr[2] / (mass_arr[1] + mass_arr[2])

    #add rotational velocity
    theta_no = rand()
    theta_no = asin(2.0 * theta_no - 1.0) + π/2.0
    vrot_no = @. [-xno[3], 0.0, xno[1]] * cos(theta_no)/sqrt(xno[1]^2 + xno[3]^2)
    vrot_no = @. vrot_no + [xno[1], -(xno[1]^2 + xno[3]^2)/xno[2], xno[3]] *
                sin(theta_no)/sqrt((xno[1]^2 + xno[3]^2) * (1 + (xno[1]^2 + xno[3]^2)/xno[2]^2))

    vvib_no = sqrt(2.0 * e_rot_i/(mass_arr[1] + mass_arr[2]))
    v0no[1, :] = @. v0no[1, :] + vvib_no * vrot_no * sqrt(mass_arr[2]/mass_arr[1])
    v0no[2, :] = @. v0no[2, :] - vvib_no * vrot_no * sqrt(mass_arr[1]/mass_arr[2])
    #stack information to x_au0 and v_au0
    v0 = [v0no; v_au0]
    x0 = [x0no; x_au0]

    #initialize nearest neigbors

    return x0, v0
end