text:
# # [Cantilever with a Tip Moment](@id tipmoment)
#
# This example shows how to predict the behavior of a cantilever beam that is subjected
# to a constant tip moment.  This is a common benchmark problem for the geometrically
# nonlinear analysis of beams.
#
# ![](../assets/tipmoment-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`tipmoment.ipynb`](@__NBVIEWER_ROOT_URL__/examples/tipmoment.ipynb).
#-

using GXBeam, LinearAlgebra

L = 12 # inches
h = w = 1 # inches
E = 30e6 # lb/in^4 Young's Modulus

A = h*w
Iyy = w*h^3/12
Izz = w^3*h/12

## bending moment (applied at end)
λ = [0.0, 0.4, 0.8, 1.2, 1.6, 1.8, 2.0]
m = pi*E*Iyy/L
M = λ*m

## create points
nelem = 16
x = range(0, L, length=nelem+1)
y = zero(x)
z = zero(x)
points = [[x[i],y[i],z[i]] for i = 1:length(x)]

## index of endpoints of each beam element
start = 1:nelem
stop = 2:nelem+1

## compliance matrix for each beam element
compliance = fill(Diagonal([1/(E*A), 0, 0, 0, 1/(E*Iyy), 1/(E*Izz)]), nelem)

## create assembly of interconnected nonlinear beams
assembly = Assembly(points, start, stop, compliance=compliance)

## pre-initialize system storage
system = System(assembly, true)

## run an analysis for each prescribed bending moment
states = Vector{AssemblyState{Float64}}(undef, length(M))
for i = 1:length(M)

    ## create dictionary of prescribed conditions
    prescribed_conditions = Dict(
        ## fixed left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        ## moment on right side
        nelem+1 => PrescribedConditions(Mz = M[i])
    )

    ## perform a static analysis
    static_analysis!(system, assembly;
        prescribed_conditions = prescribed_conditions)

    ## post-process the results
    states[i] = AssemblyState(system, assembly;
        prescribed_conditions = prescribed_conditions)

end

#!jl nothing #hide

#
# This problem has a simple analytical solution, which we obtained from "Study of the
# Geometric Stiffening Effect: Comparison of Different Formulations" by Juana M. Mayo,
# Daniel Garcia-Vallejo, and Jaime Dominguez.
#

## analytical solution (ρ = E*I/M)
analytical(x, ρ) = ifelse(ρ == Inf, zeros(3), [ρ*sin(x/ρ)-x, ρ*(1-cos(x/ρ)), 0])
#!jl nothing #hide

#
# Plotting the results reveals that the analytical and computational results show
# excellent agreement.
#

using Plots
#md using Suppressor #hide
pyplot()
#!jl nothing #hide

#-

#md @suppress_err begin #hide

## set up the plot
plot(
    xlim = (-0.25, 1.1),
    xticks = -0.25:0.25:1.0,
    xlabel = "x/L",
    ylim = (-0.05, 0.8),
    yticks = 0.0:0.2:0.8,
    ylabel = "y/L",
    aspect_ratio = 1.0,
    grid = false,
    overwrite_figure=false
    )

## create dummy legend entries for GXBeam and Analytical
scatter!([], [], color=:black, label="GXBeam")
plot!([], [], color=:black, label="Analytical")

## plot the data
for i = 1:length(M)

    local x, y

    ## GXBeam
    x = [assembly.points[ipoint][1] + states[i].points[ipoint].u[1] for ipoint =
        1:length(assembly.points)]
    y = [assembly.points[ipoint][2] + states[i].points[ipoint].u[2] for ipoint =
        1:length(assembly.points)]
    scatter!(x/L, y/L, label="", color = i)

    ## Analytical
    x0 = range(0, L, length=100)
    deflection = analytical.(x0, E*Iyy/M[i])
    x = (x0 + getindex.(deflection, 1))
    y = getindex.(deflection, 2)
    plot!(x/L, y/L, label="\$\\lambda\$=$(λ[i])", color=i)
end
#!nb plot!(show=true)
#md savefig("../assets/tipmoment-deflection.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/tipmoment-deflection.svg)

#-

# We can use this problem to test the accuracy and convergence of this package.  To do so
# we set ``\lambda = 1`` and repeat the analysis for a variety of grid sizes.  We measure
# the normalized tip displacement error ``\varepsilon(u)`` using the following expression
# ```math
# \varepsilon(u) = \left| \frac{u - u^a}{u^a} \right|
# ```
# where ``u`` is the calculated tip displacement (at x=L) and ``u^a`` is the analytical
# tip displacement.

grid_sizes = unique(round.(Int, 10 .^ range(0,3,length=25)))

L = 12 # inches
h = w = 1 # inches
E = 30e6 # lb/in^4 Young's Modulus

A = h*w
Iyy = w*h^3/12
Izz = w^3*h/12

## bending moment (applied at end)
λ = 1.0
m = pi*E*Iyy/L
M = λ*m

## run an analysis for each grid size
states = Vector{AssemblyState{Float64}}(undef, length(grid_sizes))
for (igrid, nelem) in enumerate(grid_sizes)

    local x, y, z, points, start, stop, compliance, assembly, system

    ## create points
    x = range(0, L, length=nelem+1)
    y = zero(x)
    z = zero(x)
    points = [[x[i],y[i],z[i]] for i = 1:length(x)]

    ## index of endpoints of each beam element
    start = 1:nelem
    stop = 2:nelem+1

    ## compliance matrix for each beam element
    compliance = fill(Diagonal([1/(E*A), 0, 0, 0, 1/(E*Iyy), 1/(E*Izz)]), nelem)

    ## create assembly of interconnected nonlinear beams
    assembly = Assembly(points, start, stop, compliance=compliance)

    ## create dictionary of prescribed conditions
    prescribed_conditions = Dict(
        ## fixed left side
        1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
        ## moment on right side
        nelem+1 => PrescribedConditions(Mz = M)
    )

    ## perform a static analysis
    system, converged = static_analysis(assembly;
        prescribed_conditions = prescribed_conditions)

    ## post-process the results
    states[igrid] = AssemblyState(system, assembly;
        prescribed_conditions = prescribed_conditions)

end

#!jl nothing #hide

#-

## calculate analytical solution
dxa, dya = analytical(L, E*Iyy/M)

## extract computational solution
dx = [states[igrid].points[end].u[1] for igrid = 1:length(grid_sizes)]
dy = [states[igrid].points[end].u[2] for igrid = 1:length(grid_sizes)]

## compute error
εx = abs.((dx .- dxa) ./ dxa)
εy = abs.((dy .- dya) ./ dya)

#!jl nothing #hide

#-

using Plots
#md using Suppressor #hide
pyplot()
#!jl nothing #hide

#-

#md @suppress_err begin #hide

## plot the x-error
p1 = plot(grid_sizes .+ 1, εx, label="",
    xlabel = "Number of Nodes",
    xaxis=:log,
    xlim = (10^0, 10^3),
    xtick = 10.0 .^ (0:3),
    ylabel = "\$\\varepsilon(u_x)\$",
    yaxis=:log,
    ylim = (-Inf, 10^0),
    ytick = 10.0 .^ -(0:7),
    overwrite_figure=false,
    show=true)

#!nb plot!(show=true)
#md savefig("../assets/tipmoment-x-convergence.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/tipmoment-x-convergence.svg)

#-

#md @suppress_err begin #hide

## plot the y-error
p2 = plot(grid_sizes .+ 1, εy, label="",
    xlabel = "Number of Nodes",
    xaxis=:log,
    xlim = (10^0, 10^3),
    xtick = 10.0 .^ (0:3),
    ylabel = "\$\\varepsilon(u_y)\$",
    yaxis=:log,
    ylim = (-Inf, 10^0),
    ytick = 10.0 .^ -(0:7),
    overwrite_figure=false,
    show=true)

#!nb plot!(show=true)
#md savefig("../assets/tipmoment-y-convergence.svg") #hide
#md closeall() #hide
#md end #hide
#md nothing #hide

#md # ![](../assets/tipmoment-y-convergence.svg)

#-

# We observe second-order algebraic convergence for both x and y tip displacement errors.
# We can therefore conclude that a large number of elements are likely necessary in order
# to obtain highly accurate solutions using this package.  For problems where high
# accuracy solutions are critical, higher order shape functions, such as the Legendre
# spectral finite elements used by [BeamDyn](https://www.nrel.gov/wind/nwtc/beamdyn.html)
# are likely more computationally efficient.


text:
#A counter type

struct Accumulator{T, V<:Number} <: AbstractDict{T,V}
    map::Dict{T,V}
end

## constructors

Accumulator{T, V}() where {T,V<:Number} = Accumulator{T,V}(Dict{T,V}())
@deprecate Accumulator(::Type{T}, ::Type{V}) where {T,V<:Number} Accumulator{T, V}()

counter(T::Type) = Accumulator{T,Int}()

counter(dct::Dict{T,V}) where {T,V<:Integer} = Accumulator{T,V}(copy(dct))

"""
    counter(seq)

Returns an `Accumulator` object containing the elements from `seq`.
"""
function counter(seq)
    ct = counter(eltype_for_accumulator(seq))
    for x in seq
        inc!(ct, x)
    end
    return ct
end

eltype_for_accumulator(seq::T) where T = eltype(T)
function eltype_for_accumulator(seq::T) where {T<:Base.Generator}
    Base.@default_eltype(seq)
end

copy(ct::Accumulator) = Accumulator(copy(ct.map))

length(a::Accumulator) = length(a.map)

## retrieval

get(ct::Accumulator, x, default) = get(ct.map, x, default)
# need to allow user specified default in order to
# correctly implement "informal" AbstractDict interface

getindex(ct::Accumulator{T,V}, x) where {T,V} = get(ct.map, x, zero(V))

setindex!(ct::Accumulator, x, v) = setindex!(ct.map, x, v)

haskey(ct::Accumulator, x) = haskey(ct.map, x)

keys(ct::Accumulator) = keys(ct.map)

values(ct::Accumulator) = values(ct.map)

sum(ct::Accumulator) = sum(values(ct.map))

## iteration

iterate(ct::Accumulator, s...) = iterate(ct.map, s...)

# manipulation

"""
    inc!(ct, x, [v=1])

Increments the count for `x` by `v` (defaulting to one)
"""
inc!(ct::Accumulator, x, a::Number) = (ct[x] += a)
inc!(ct::Accumulator{T,V}, x) where {T,V} = inc!(ct, x, one(V))

# inc! is preferred over push!, but we need to provide push! for the Bag interpreation
# which is used by classified_collections.jl
push!(ct::Accumulator, x) = inc!(ct, x)
push!(ct::Accumulator, x, a::Number) = inc!(ct, x, a)

# To remove ambiguities related to Accumulator now being a subtype of AbstractDict
push!(ct::Accumulator, x::Pair)  = inc!(ct, x)

"""
    dec!(ct, x, [v=1])

Decrements the count for `x` by `v` (defaulting to one)
"""
dec!(ct::Accumulator, x, a::Number) = (ct[x] -= a)
dec!(ct::Accumulator{T,V}, x) where {T,V} = dec!(ct, x, one(V))

#TODO: once we are done deprecating `pop!` for `reset!` then add `pop!` as an alias for `dec!`

"""
    merge!(ct1, others...)

Merges the other counters into `ctl`,
summing the counts for all elements.
"""
function merge!(ct::Accumulator, other::Accumulator)
    for (x, v) in other
        inc!(ct, x, v)
    end
    ct
end

function merge!(ct1::Accumulator, others::Accumulator...)
    for ct in others
        merge!(ct1,ct)
    end
    return ct1
end

"""
     merge(counters...)

Creates a new counter with total counts equal to the sum of the counts in the counters given as arguments.

See also merge!
"""
function merge(ct1::Accumulator, others::Accumulator...)
    ct = copy(ct1)
    merge!(ct,others...)
end

"""
    reset!(ct::Accumulator, x)

Resets the count of `x` to zero.
Returns its former count.
"""
reset!(ct::Accumulator, x) = pop!(ct.map, x)

"""
     nlargest(acc::Accumulator, [n])

Returns a sorted vector of the `n` most common elements, with their counts.
If `n` is omitted, the full sorted collection is returned.

This corresponds to Python's `Counter.most_common` function.

Example
```
julia> nlargest(counter("abbbccddddda"))

4-element Array{Pair{Char,Int64},1}:
 'd'=>5
 'b'=>3
 'c'=>2
 'a'=>2

julia> nlargest(counter("abbbccddddda"),2)

2-element Array{Pair{Char,Int64},1}:
 'd'=>5
 'b'=>3

```
"""
nlargest(acc::Accumulator) = sort!(collect(acc), by=last, rev=true)
nlargest(acc::Accumulator, n) = partialsort!(collect(acc), 1:n, by=last, rev=true)

"""
     nsmallest(acc::Accumulator, [n])

Returns a sorted vector of the `n` least common elements, with their counts.
If `n` is omitted, the full sorted collection is returned.

This is the opposite of the `nlargest` function.
For obvious reasons this will not include zero counts for items not encountered.
(unless those elements are added to he accumulator directly, eg via `acc[foo]=0)
"""
nsmallest(acc::Accumulator) = sort!(collect(acc), by=last, rev=false)
nsmallest(acc::Accumulator, n) = partialsort!(collect(acc), 1:n, by=last, rev=false)

## Deprecations
@deprecate pop!(ct::Accumulator, x) reset!(ct, x)
@deprecate push!(ct1::Accumulator, ct2::Accumulator) merge!(ct1,ct2)

ext: jl


text:
using Printf

function output_restart(stepnum,nodes,Qbase,cellxmax,cellymax,out_file_front,out_ext,out_dir)

    stepnum=string(stepnum)
    while length(stepnum) < 6
        stepnum="0"*stepnum
    end

    fff=out_dir*"/"*out_file_front*stepnum*out_ext
    open(fff,"w") do f
        write(f,"result:x,y,T[K]\n")
        for i in 2:cellxmax-1
            for j in 2:cellymax-1
                # x = 0.25*(nodes[i+1,j+1,1]+nodes[i+1,j,1]+nodes[i,j+1,1]+nodes[i,j,1])
                # y = 0.25*(nodes[i+1,j+1,2]+nodes[i+1,j,2]+nodes[i,j+1,2]+nodes[i,j,2])
                # write(f,string(x)*" "*string(y)*" "*string(Qbase[i,j])*"\n")
                a = @sprintf("%8.8f", Qbase[i,j])
                write(f, a*"\n")
            end
        end
    end
    println("write "*fff)
end

function output_result(stepnum,nodes,Qbase,cellxmax,cellymax,out_file_front,out_ext,out_dir)

    stepnum=string(stepnum)
    while length(stepnum) < 6
        stepnum="0"*stepnum
    end

    fff=out_dir*"/"*out_file_front*stepnum*out_ext
    open(fff,"w") do f
        write(f,"result:x,y,T[K]\n")
        for i in 2:cellxmax-1
            for j in 2:cellymax-1
                # x = 0.25*(nodes[i+1,j+1,1]+nodes[i+1,j,1]+nodes[i,j+1,1]+nodes[i,j,1])
                # y = 0.25*(nodes[i+1,j+1,2]+nodes[i+1,j,2]+nodes[i,j+1,2]+nodes[i,j,2])
                # write(f,string(x)*" "*string(y)*" "*string(Qbase[i,j])*"\n")
                a = @sprintf("%8.8f", Qbase[i,j])
                write(f, a*"\n")
            end
        end
    end
    println("\nwrite "*fff)
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct CSVInput <: File
        name::String
        descriptor::String
    end

Read .csv file

# Arguments
- `name::String`: input file name
- `descriptor::String`: file descriptor
"""
mutable struct CSVInput <: File
    "input file name"
    name::String
    "file descriptor"
    descriptor::String
end

function CSVInput(; name, descriptor, )
    CSVInput(name, descriptor, )
end

"""Get [`CSVInput`](@ref) `name`."""
get_name(value::CSVInput) = value.name
"""Get [`CSVInput`](@ref) `descriptor`."""
get_descriptor(value::CSVInput) = value.descriptor

"""Set [`CSVInput`](@ref) `name`."""
set_name!(value::CSVInput, val) = value.name = val
"""Set [`CSVInput`](@ref) `descriptor`."""
set_descriptor!(value::CSVInput, val) = value.descriptor = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using Gadfly, DataArrays, RDatasets

set_default_plot_size(6inch, 3inch)

plot(dataset("ggplot2", "diamonds"),
     x=:Price, color=:Cut,
     Geom.histogram)


text:
mktempdir() do tmpdir
    # Test dataset multiplication
    n = 30
    m = 10
    a = rand(n,n) + im * rand(n,n)

    fn = joinpath(tmpdir, "exactdiag-mul.jld")
    JLD.save(fn, "a", a)

    jldopen(fn) do file
        b = rand(n, m) + im * rand(n, m)
        @test ExactDiag.my_Ac_mul_B(file["a"], b) ≈ a' * b

        c = Array{ComplexF64}(undef, size(b)...)
        @test ExactDiag.my_A_mul_B!(c, file["a"], b) ≈ a * b
    end
end

let
    lattice = ChainLattice([8])
    apply_hamiltonian = spin_half_hamiltonian(J1=1)
    hs = SpinHalfHilbertSpace(lattice)
    seed_state!(hs, N_up=div(length(lattice), 2))
    rst = RepresentativeStateTable(hs, apply_hamiltonian)

    initial_state = zeros(ComplexF64, length(rst.hs.indexer))
    initial_state[1] = 1 / sqrt(2)
    initial_state[2] = im / sqrt(2)
    time_steps = 10 .^ range(-1.5, stop=4, length=41)
    push!(time_steps, 0)

    function calculate_momentum_sector(func, sector_index, momentum_index)
        diagsect = DiagonalizationSector(rst, sector_index, momentum_index)
        reduced_hamiltonian = Matrix(construct_reduced_hamiltonian(diagsect))
        fact = eigen(Hermitian((reduced_hamiltonian + reduced_hamiltonian') / 2))
        func(construct_reduced_indexer(diagsect), fact.values, fact.vectors)
        nothing
    end

    # Test that the initial state can be either a Vector or Matrix, with the same results
    let eb = ExactDiag.to_energy_basis(calculate_momentum_sector, rst, initial_state)[1]
        @test ExactDiag.to_energy_basis(calculate_momentum_sector, rst, [initial_state initial_state])[1] == [eb eb]
    end

    let output_states = time_evolve(calculate_momentum_sector, rst, initial_state, time_steps)
        # Test that time-evolved norms are all (essentially) unity
        for t_i in 1:length(time_steps)
            @test norm(output_states[:, t_i]) ≈ 1
        end

        # Test that evolving for zero time returns that same wavefunction
        @test initial_state ≈ output_states[:, end]

        # Test that evolving forward, then backwards, returns the same wavefunction
        reverse_evolved_state = time_evolve(calculate_momentum_sector, rst, output_states[:, end], [-time_steps[end]])
        @test initial_state ≈ reverse_evolved_state

        # Test JLD dataset time evolution
        mktempdir() do tmpdir
            fn = joinpath(tmpdir, "exactdiag-evolve.jld")
            jldopen(fn, "w") do file
                # Save each momentum sector
                for momentum_index in eachmomentumindex(lattice)
                    calculate_momentum_sector(1, momentum_index) do indexer, evals, evecs
                        file["indexer_$(momentum_index)"] = indexer
                        file["evals_$(momentum_index)"] = evals
                        file["evecs_$(momentum_index)"] = evecs
                        nothing
                    end
                end

                # Perform time evolution using the saved `JldDataset`s
                output_states_jld = time_evolve(rst, initial_state, time_steps) do func, sector_index, momentum_index
                    @test sector_index == 1
                    indexer = read(file["indexer_$(momentum_index)"])
                    evals = read(file["evals_$(momentum_index)"])
                    evecs = file["evecs_$(momentum_index)"]
                    func(indexer, evals, evecs)
                    nothing
                end

                # Check the results
                @test output_states ≈ output_states_jld
            end
        end

        # Test multiple initial states
        let output_states2 = time_evolve(calculate_momentum_sector, rst, [initial_state (im * initial_state)], time_steps)
            @test output_states2[:, :, 1] ≈ output_states
            @test output_states2[:, :, 2] ≈ im * output_states
        end

        # Test that doing exact diagonalization without regard to momentum sectors
        # would give the same results
        full_ham = operator_matrix(hs, apply_hamiltonian)
        fact = eigen(Hermitian(Matrix(full_ham)))
        ψ_e = fact.vectors' * initial_state
        ψ_t = fact.vectors * (exp.(-im .* fact.values .* transpose(time_steps)) .* ψ_e)
        @test ψ_t ≈ output_states
    end

    # Try evolving a state whose support is only on a subset of momentum sectors
    initial_state = zeros(ComplexF64, length(rst.hs.indexer))
    initial_state[findfirst(isequal([0,1,0,1,0,1,0,1]), hs.indexer)] = 1

    let output_states = time_evolve(calculate_momentum_sector, rst, initial_state, time_steps, k_indices=[1,5])
        # Test the output is of the correct size
        @test size(output_states) == (length(rst.hs.indexer), length(time_steps))

        # Test that time-evolved norms are all (essentially) unity
        for t_i in 1:length(time_steps)
            @test norm(output_states[:, t_i]) ≈ 1
        end

        # Test that evolving for zero time returns that same wavefunction
        @test initial_state ≈ output_states[:, end]

        # Test against standard evolution
        let output_states2 = time_evolve(calculate_momentum_sector, rst, initial_state, time_steps)
            @test output_states2 ≈ output_states
        end
    end

    # Try providing an incomplete subset of momentum sectors
    let
        ψ_e, = ExactDiag.to_energy_basis(calculate_momentum_sector, rst, initial_state, k_indices=[1])

        # Test that the vector is the correct size
        @test length(ψ_e) < length(rst.hs.indexer)

        # Test that the norm is not unity.
        @test norm(ψ_e) ≈ 1/sqrt(2)
    end
end


text:
# Julia wrapper for header: /usr/include/libavutil/channel_layout.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_get_channel_layout,
    av_get_extended_channel_layout,
    av_get_channel_layout_string,
    av_bprint_channel_layout,
    av_get_channel_layout_nb_channels,
    av_get_default_channel_layout,
    av_get_channel_layout_channel_index,
    av_channel_layout_extract_channel,
    av_get_channel_name,
    av_get_channel_description,
    av_get_standard_channel_layout

function av_get_channel_layout(name)
    ccall((:av_get_channel_layout, libavutil), UInt64, (Cstring,), name)
end

function av_get_extended_channel_layout(name, channel_layout, nb_channels)
    ccall((:av_get_extended_channel_layout, libavutil), Cint, (Cstring, Ptr{UInt64}, Ptr{Cint}), name, channel_layout, nb_channels)
end

function av_get_channel_layout_string(buf, buf_size::Integer, nb_channels::Integer, channel_layout::UInt64)
    ccall((:av_get_channel_layout_string, libavutil), Cvoid, (Cstring, Cint, Cint, UInt64), buf, buf_size, nb_channels, channel_layout)
end

function av_bprint_channel_layout(bp, nb_channels::Integer, channel_layout::UInt64)
    ccall((:av_bprint_channel_layout, libavutil), Cvoid, (Ptr{AVBPrint}, Cint, UInt64), bp, nb_channels, channel_layout)
end

function av_get_channel_layout_nb_channels(channel_layout::UInt64)
    ccall((:av_get_channel_layout_nb_channels, libavutil), Cint, (UInt64,), channel_layout)
end

function av_get_default_channel_layout(nb_channels::Integer)
    ccall((:av_get_default_channel_layout, libavutil), Int64, (Cint,), nb_channels)
end

function av_get_channel_layout_channel_index(channel_layout::UInt64, channel::UInt64)
    ccall((:av_get_channel_layout_channel_index, libavutil), Cint, (UInt64, UInt64), channel_layout, channel)
end

function av_channel_layout_extract_channel(channel_layout::UInt64, index::Integer)
    ccall((:av_channel_layout_extract_channel, libavutil), UInt64, (UInt64, Cint), channel_layout, index)
end

function av_get_channel_name(channel::UInt64)
    ccall((:av_get_channel_name, libavutil), Cstring, (UInt64,), channel)
end

function av_get_channel_description(channel::UInt64)
    ccall((:av_get_channel_description, libavutil), Cstring, (UInt64,), channel)
end

function av_get_standard_channel_layout(index::Integer, layout, name)
    ccall((:av_get_standard_channel_layout, libavutil), Cint, (UInt32, Ptr{UInt64}, Ptr{Cstring}), index, layout, name)
end


text:
"""
	lsa(dtm::DocumentTermMatrix)
	lsa(crps::Corpus)

Performs Latent Semantic Analysis or LSA on a corpus.

"""
lsa(dtm::DocumentTermMatrix) = svd(Matrix(tf_idf(dtm)))
function lsa(crps::Corpus)
	update_lexicon!(crps)
	svd(Matrix(tf_idf(DocumentTermMatrix(crps))))
end


text:
# ---
# title: 983. Minimum Cost For Tickets
# id: problem983
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/minimum-cost-for-tickets/description/>
# hidden: true
# ---
#
# In a country popular for train travel, you have planned some train travelling
# one year in advance.  The days of the year that you will travel is given as an
# array `days`.  Each day is an integer from `1` to `365`.
#
# Train tickets are sold in 3 different ways:
#
#   * a 1-day pass is sold for `costs[0]` dollars;
#   * a 7-day pass is sold for `costs[1]` dollars;
#   * a 30-day pass is sold for `costs[2]` dollars.
#
# The passes allow that many days of consecutive travel.  For example, if we get
# a 7-day pass on day 2, then we can travel for 7 days: day 2, 3, 4, 5, 6, 7,
# and 8.
#
# Return the minimum number of dollars you need to travel every day in the given
# list of `days`.
#
#
#
# **Example 1:**
#
#
#
#     Input: days = [1,4,6,7,8,20], costs = [2,7,15]
#     Output: 11
#     Explanation:
#     For example, here is one way to buy passes that lets you travel your travel plan:
#     On day 1, you bought a 1-day pass for costs[0] = /$/2, which covered day 1.
#     On day 3, you bought a 7-day pass for costs[1] = /$/7, which covered days 3, 4, ..., 9.
#     On day 20, you bought a 1-day pass for costs[0] = /$/2, which covered day 20.
#     In total you spent /$/11 and covered all the days of your travel.
#
#
# **Example 2:**
#
#
#
#     Input: days = [1,2,3,4,5,6,7,8,9,10,30,31], costs = [2,7,15]
#     Output: 17
#     Explanation:
#     For example, here is one way to buy passes that lets you travel your travel plan:
#     On day 1, you bought a 30-day pass for costs[2] = /$/15 which covered days 1, 2, ..., 30.
#     On day 31, you bought a 1-day pass for costs[0] = /$/2 which covered day 31.
#     In total you spent /$/17 and covered all the days of your travel.
#
#
#
#
# **Note:**
#
#   1. `1 <= days.length <= 365`
#   2. `1 <= days[i] <= 365`
#   3. `days` is in strictly increasing order.
#   4. `costs.length == 3`
#   5. `1 <= costs[i] <= 1000`
#
#
## @lc code=start
using LeetCode

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


text:
module Cosmology

export cosmology,
       age_gyr,
       angular_diameter_dist_mpc,
       comoving_radial_dist_mpc,
       comoving_transverse_dist_mpc,
       comoving_volume_gpc3,
       H,
       hubble_dist_mpc,
       hubble_time_gyr,
       luminosity_dist_mpc,
       lookback_time_gyr,
       scale_factor

abstract AbstractCosmology
abstract AbstractClosedCosmology <: AbstractCosmology
abstract AbstractFlatCosmology <: AbstractCosmology
abstract AbstractOpenCosmology <: AbstractCosmology

immutable FlatLCDM <: AbstractFlatCosmology
    h::Float64
    Ω_Λ::Float64
    Ω_m::Float64
    Ω_r::Float64

    function FlatLCDM(h::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real)
        new(float64(h), float64(Ω_Λ), float64(Ω_m), float64(Ω_r))
    end
end

a2E(c::FlatLCDM, a::Float64) = sqrt(c.Ω_r + c.Ω_m*a + c.Ω_Λ*a^4)

immutable ClosedLCDM <: AbstractClosedCosmology
    h::Float64
    Ω_k::Float64
    Ω_Λ::Float64
    Ω_m::Float64
    Ω_r::Float64

    function ClosedLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real)
        new(float64(h), float64(Ω_k), float64(Ω_Λ), float64(Ω_m), float64(Ω_r))
    end
end

immutable OpenLCDM <: AbstractOpenCosmology
    h::Float64
    Ω_k::Float64
    Ω_Λ::Float64
    Ω_m::Float64
    Ω_r::Float64

    function OpenLCDM(h::Real, Ω_k::Real, Ω_Λ::Real, Ω_m::Real, Ω_r::Real)
        new(float64(h), float64(Ω_k), float64(Ω_Λ), float64(Ω_m), float64(Ω_r))
    end
end

function a2E(c::Union(ClosedLCDM,OpenLCDM), a::Float64)
    a2 = a*a
    sqrt(c.Ω_r + c.Ω_m*a + (c.Ω_k + c.Ω_Λ*a2)*a2)
end

for c in ("Flat", "Open", "Closed")
    @eval begin
        immutable $(symbol("$(c)WCDM")) <: $(symbol("Abstract$(c)Cosmology"))
            h::Float64
            Ω_k::Float64
            Ω_Λ::Float64
            Ω_m::Float64
            Ω_r::Float64
            w0::Float64
            wa::Float64
        end
    end
end

function WCDM(h::Float64, Ω_k::Float64, Ω_Λ::Float64, Ω_m::Float64, Ω_r::Float64, w0::Float64, wa::Float64)
    if Ω_k < 0
        ClosedWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    elseif Ω_k > 0
        OpenWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    else
        FlatWCDM(h, Ω_k, Ω_Λ, Ω_m, Ω_r, w0, wa)
    end
end

function a2E(c::Union(FlatWCDM,ClosedWCDM,OpenWCDM), a::Float64)
    ade = exp((1.0 - 3.0*(c.w0 + c.wa))*log(a) + 3.0*c.wa*(a - 1.0))
    sqrt(c.Ω_r + (c.Ω_m + c.Ω_k*a)*a + c.Ω_Λ*ade)
end

function cosmology(;h=0.69,
                   Neff=3.04,
                   OmegaK=0,
                   OmegaM=0.29,
                   OmegaR=nothing,
                   Tcmb=2.7255,
                   w0=-1,
                   wa=0)

    if OmegaR === nothing
        OmegaG = 4.48131e-7*Tcmb^4/h^2
        OmegaN = Neff*OmegaG*(7/8)*(4/11)^(4/3)
        OmegaR = OmegaG + OmegaN
    end

    OmegaL = 1. - OmegaK - OmegaM - OmegaR

    if !(w0 == -1 && wa == 0)
        return WCDM(float64(h), float64(OmegaK), float64(OmegaL),
                    float64(OmegaM), float64(OmegaR), float64(w0), float64(wa))
    end

    if OmegaK < 0
        return ClosedLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
    elseif OmegaK > 0
        return OpenLCDM(h, OmegaK, OmegaL, OmegaM, OmegaR)
    else
        return FlatLCDM(h, OmegaL, OmegaM, OmegaR)
    end
end

# hubble rate

scale_factor(z) = 1/(1 + z)
E(c::AbstractCosmology, z) = (a = scale_factor(z); a2E(c,a)/a^2)
H(c::AbstractCosmology, z) = 100. * c.h * E(c, z)

hubble_dist_mpc0(c::AbstractCosmology) = 2997.92458/c.h
hubble_dist_mpc(c::AbstractCosmology, z) = hubble_dist_mpc0(c)/E(c,z)

hubble_time_gyr0(c::AbstractCosmology) = 9.77814/c.h
hubble_time_gyr(c::AbstractCosmology, z) = hubble_time_gyr0(c)/E(c,z)

# distances

Z(c::AbstractCosmology, z::Real) = ((q,_) = quadgk(a::Float64->1.0/a2E(c,a), scale_factor(z), 1); q)

comoving_radial_dist_mpc(c::AbstractCosmology, z) = hubble_dist_mpc0(c)*Z(c, z)

comoving_transverse_dist_mpc(c::AbstractFlatCosmology, z) =
    comoving_radial_dist_mpc(c, z)
function comoving_transverse_dist_mpc(c::AbstractOpenCosmology, z)
    sqrtok = sqrt(c.Ω_k)
    hubble_dist_mpc0(c)*sinh(sqrtok*Z(c,z))/sqrtok
end
function comoving_transverse_dist_mpc(c::AbstractClosedCosmology, z)
    sqrtok = sqrt(abs(c.Ω_k))
    hubble_dist_mpc0(c)*sin(sqrtok*Z(c,z))/sqrtok
end

angular_diameter_dist_mpc(c::AbstractCosmology, z) =
    comoving_transverse_dist_mpc(c, z)/(1 + z)

luminosity_dist_mpc(c::AbstractCosmology, z) =
    comoving_transverse_dist_mpc(c, z)*(1 + z)

# volumes

comoving_volume_gpc3(c::AbstractFlatCosmology, z) =
    (4pi/3)*(comoving_radial_dist_mpc(c,z)*1e-3)^3
function comoving_volume_gpc3(c::AbstractOpenCosmology, z)
    DH = hubble_dist_mpc0(c)
    x = comoving_transverse_dist_mpc(c,z)/DH
    sqrtok = sqrt(c.Ω_k)
    2pi*(DH*1e-3)^3*(x*sqrt(1. + c.Ω_k*x^2) - asinh(sqrtok*x)/sqrtok)/c.Ω_k
end
function comoving_volume_gpc3(c::AbstractClosedCosmology, z)
    DH = hubble_dist_mpc0(c)
    x = comoving_transverse_dist_mpc(c,z)/DH
    sqrtok = sqrt(abs(c.Ω_k))
    2pi*(DH*1e-3)^3*(x*sqrt(1. + c.Ω_k*x^2) - asin(sqrtok*x)/sqrtok)/c.Ω_k
end

comoving_volume_element_gpc3(c::AbstractCosmology, z) =
    1e-9*hubble_dist_mpc0(c,z)*angular_diameter_dist_mpc(c,z)^2/a2E(c,scale_factor(z))

# times

T(c::AbstractCosmology, a0::Float64, a1::Float64) = ((q,_) = quadgk(x::Float64->x/a2E(c,x), a0, a1); q)
age_gyr(c::AbstractCosmology, z) = hubble_time_gyr0(c)*T(c, 0., scale_factor(z))
lookback_time_gyr(c::AbstractCosmology, z) = hubble_time_gyr0(c)*T(c, scale_factor(z), 1.)

end # module


text:
"bounding box representing the smallest straightup rectangle containing all data "
struct BoundingBox
    xmin::Float64
    ymin::Float64
    zmin::Float64
    xmax::Float64
    ymax::Float64
    zmax::Float64
end

# construct BoundingBox with keyword argumets to prevent mixups
function BoundingBox(;xmin=0.0::Float64, ymin=0.0::Float64, zmin=0.0::Float64,
    xmax=0.0::Float64, ymax=0.0::Float64, zmax=0.0::Float64)
    BoundingBox(xmin, ymin, zmin, xmax, ymax, zmax)
end

function Base.show(io::IO, bbox::BoundingBox)
    print(io, @sprintf "xmin=%.1f, " bbox.xmin )
    print(io, @sprintf "xmax=%.1f, " bbox.xmax )
    print(io, @sprintf "ymin=%.1f, " bbox.ymin )
    print(io, @sprintf "ymax=%.1f, " bbox.ymax )
    print(io, @sprintf "zmin=%.2f, " bbox.zmin )
    print(io, @sprintf "zmax=%.2f" bbox.zmax )
end

"""Cloud is a point cloud container type similar to Clouds' Cloud, without
a spatial index, but with a simple bounding box that is kept up to date."""
mutable struct Cloud
    positions::Vector{MVector{3,Float64}}
    attributes::Dict{Symbol,Vector}
    boundingbox::BoundingBox
end

Cloud(pos, attrs) = Cloud(pos, attrs, calc_bbox(pos))

function Base.show(io::IO, cloud::Cloud)
    print(io, "PointCloud (N=$(length(cloud)), bbox=[$(cloud.boundingbox)], attributes=$(keys(cloud.attributes)))")
end

"Get the vector of point cloud positions"
positions(cloud::Cloud) = getfield(cloud, :positions)

"Get the Dict of point cloud attributes"
attributes(cloud::Cloud) = getfield(cloud, :attributes)

"Get the bounding box of the point cloud"
boundingbox(cloud::Cloud) = getfield(cloud, :boundingbox)

"Calculate bounding box of the positions"
function calc_bbox(P::Vector{SVector{3,Float64}})
    xmin, xmax, ymin, ymax, zmin, zmax = Inf, -Inf, Inf, -Inf, Inf, -Inf
    @inbounds for p in P
        xmin = min(xmin, p[1])
        xmax = max(xmax, p[1])
        ymin = min(ymin, p[2])
        ymax = max(ymax, p[2])
        zmin = min(zmin, p[3])
        zmax = max(zmax, p[3])
    end
    BoundingBox(xmin=xmin,ymin=ymin,zmin=zmin,xmax=xmax,ymax=ymax,zmax=zmax)
end

"Calculate bounding box of the positions"
function calc_bbox(P::Vector{MVector{3,Float64}})
    xmin, xmax, ymin, ymax, zmin, zmax = Inf, -Inf, Inf, -Inf, Inf, -Inf
    @inbounds for p in P
        xmin = min(xmin, p[1])
        xmax = max(xmax, p[1])
        ymin = min(ymin, p[2])
        ymax = max(ymax, p[2])
        zmin = min(zmin, p[3])
        zmax = max(zmax, p[3])
    end
    BoundingBox(xmin=xmin,ymin=ymin,zmin=zmin,xmax=xmax,ymax=ymax,zmax=zmax)
end

"Calculate bounding box of the point cloud"
function calc_bbox(cloud::Cloud)
    P = positions(cloud)
    calc_bbox(P)
end

"Update the bounding box of the point cloud"
function update_bbox!(cloud::Cloud)
    cloud.boundingbox = calc_bbox(cloud)
    nothing
end

"Subset a point cloud with specified indices"
function Base.getindex(cloud::Cloud, I)
    pos = positions(cloud)[I]
    attrs = Dict{Symbol,Vector}()
    for (k,v) in attributes(cloud)
        attrs[k] = v[I]
    end
    Cloud(pos, attrs)
end

"Get a point cloud attribute by indexing with a `Symbol`"
Base.getindex(cloud::Cloud, attr::Symbol) = attributes(cloud)[attr]

"Set a point cloud attribute assigning to a `Symbol` index"
Base.setindex!(cloud::Cloud, A::AbstractVector, attr::Symbol) = attributes(cloud)[attr] = A

"Number of points in the point cloud"
Base.length(cloud::Cloud) = length(positions(cloud))


text:
using Documenter

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

using StippleUI

makedocs(
    sitename = "SearchLight - Concise, secure, cross-platform query builder and ORM for Julia",
    format = Documenter.HTML(prettyurls = false),
    pages = [
        "Home" => "index.md",
        "StippleUI API" => [
          "API" => "api/API.md",
          "Badges" => "api/badges.md",
          "Banners" => "api/banners.md",
          "BigNumbers" => "api/bignumbers.md",
          "Buttons" => "api/buttons.md",
          "Cards" => "api/cards.md",
          "Checkboxes" => "api/checkboxes.md",
          "Chips" => "api/chips.md",
          "Dashboards" => "api/dashboards.md",
          "DatePickers" => "api/datepickers.md",
          "Dialogs" => "api/dialogs.md",
          "Drawers" => "api/drawers.md",
          "Editors" => "api/editors.md",
          "FormInputs" => "api/forminputs.md",
          "Forms" => "api/forms.md",
          "Headings" => "api/headings.md",
          "Icons" => "api/icons.md",
          "InnerLoaders" => "api/innerloaders.md",
          "Intersections" => "api/intersections.md",
          "Knobs" => "api/knobs.md",
          "Layouts" => "api/layouts.md",
          "Lists" => "api/lists.md",
          "Menus" => "api/menus.md",
          "PopupProxies" => "api/popupproxies.md",
          "Radios" => "api/radios.md",
          "Ranges" => "api/ranges.md",
          "ScrollAreas" => "api/scrollareas.md",
          "Selects" => "api/selects.md",
          "Separators" => "api/separators.md",
          "Spaces" => "api/spaces.md",
          "Spinners" => "api/spinners.md",
          "StippleUI" => "api/stippleui.md",
          "Tables" => "api/tables.md",
          "Toggles" => "api/toggles.md",
          "Uploaders" => "api/uploaders.md",
        ]
    ],
)

deploydocs(
  repo = "github.com/GenieFramework/StippleUI.jl.git",
)


text:
using EngEconomics

# Given
aInit = 2000
aLife = 2
aSalvage = 0
aAnnualSavings = 1200
bInit = 4000
bLife = 3
bSalvage = 3000
bAnnualSavingA = 2000
bAnnualSavingsG = -200
cInit = 5000
cLife = 6
cSalvage = 1000
cAnnualSavingsA = 1400
cAnnualSavingsg = 0.05
i = 0.07

# Find
years = 6
aPeriods = years / aLife
bPeriods = years / bLife
cPeriods = years / cLife

# Device A
aInitPW = -aInit * (1 + presentWorthFactor(i, aLife) + presentWorthFactor(i, 2 * aLife))
aAnnualBenefitPW = (aAnnualBenefit * seriesPresentAmountFactor(i, aLife)) * (1 + presentWorthFactor(i, aLife) + presentWorthFactor(i, 2 * aLife))
aSalvagePW = 0
aPW = aInitPW + aAnnualBenefitPW + aSalvagePW

bInitPW = -bInit * (1 + presentWorthFactor(i, bLife))
bAnnualBenefitPW = (bAnnualSavings * seriesPresentAmountFactor(i, bLife) + bAnnualSavingsG * singlePaymentPresentWorthArithmetic(i, bLife)) * (1 + presentWorthFactor(i, bLife))
bSalvagePW = bSalvage * (presentWorthFactor(i, bLife) + presentWorthFactor(i, 2 * bLife))
bPW = bInitPW + bAnnualBenefitPW + bSalvagePW

cInitPW = -cInit
cAnnualBenefitPW = (cAnnualSavingsA * presentWorthConversionFactorGeometric(cAnnualSavingsg, i, cLife))
cSalvagePW = cSalvage * presentWorthFactor(i, cLife)
cPW = cInitPW + cAnnualBenefitPW + cSalvagePW

newBSalvage = (cPW - bInitPW + bAnnualBenefitPW) / (presentWorthFactor(i, bLife) + presentWorthFactor(i, 2 * bLife))


text:
include(joinpath("src", "argparse.jl"))
const N, K, upper_bound, scseps = let
    args = parse_commandline()
    args["n"], args["k"], args["lambda"], args["tol"]
end

using SparseArrays
using LinearAlgebra
using Dates

using AbstractAlgebra
using PropertyT
using PropertyT.Groups
using PropertyT.GroupRings
using PropertyT.IntervalArithmetic

using SCS
import PropertyT.JuMP
using PropertyT.JLD

@info "Running checks for Adj_$N + $K·Op_$N - $upper_bound·Δ_$N"

G = AutGroup(FreeGroup(N), special=true)
S = PropertyT.generating_set(G)

const prefix = "SAutF$(N)_r2"
isdir(prefix) || mkpath(prefix)

const DELTA_FILE = joinpath(prefix,"delta.jld")
const SQADJOP_FILE = joinpath(prefix, "SqAdjOp_coeffs.jld")
const BLOCKDEC_FILE = joinpath(prefix, "BlockDecomposition.jld")

const fullpath = joinpath(prefix, "$(upper_bound)_K=$K")
isdir(fullpath) || mkpath(fullpath)
const WARMSTART_FILE = joinpath(fullpath, "warmstart.jld")
const SOLUTION_FILE = joinpath(fullpath, "solution.jld")

@info "Looking for delta.jld, SqAdjOp_coeffs.jld and BlockDecomposition.jld in $prefix"

if isfile(DELTA_FILE) && isfile(SQADJOP_FILE) && isfile(BLOCKDEC_FILE)
    # cached
    Δ = PropertyT.loadGRElem(DELTA_FILE, G)
    RG = parent(Δ)
    block_decomposition = load(BLOCKDEC_FILE, "BlockDecomposition")
    sq_c, adj_c, op_c = load(SQADJOP_FILE, "Sq", "Adj", "Op")
    sq = GroupRingElem(sq_c, RG)
    adj = GroupRingElem(adj_c, RG)
    op = GroupRingElem(op_c, RG);
else
    @info "Computing Laplacian"
    Δ = PropertyT.Laplacian(S, 2)
    PropertyT.saveGRElem(DELTA_FILE, Δ)
    RG = parent(Δ)

    @info "Computing Sq, Adj, Op"
    @time sq, adj, op = PropertyT.SqAdjOp(RG, N)

    save(SQADJOP_FILE, "Sq", sq.coeffs, "Adj", adj.coeffs, "Op", op.coeffs)

    @info "Compute BlockDecomposition"
    if !isfile(BLOCKDEC_FILE)
        block_decomposition = PropertyT.BlockDecomposition(RG, WreathProduct(SymmetricGroup(2), SymmetricGroup(N)))
        save(BLOCKDEC_FILE, "BlockDecomposition", block_decomposition)
    else
        block_decomposition = load(BLOCKDEC_FILE, "BlockDecomposition")
    end
end;

block_decomposition = PropertyT.decimate(block_decomposition);

elt = adj + K*op;
ELT_STRING = "Adj_$(N)+$(K)·Op_$(N)"

@info "Looking for solution.jld in $fullpath"

if !isfile(SOLUTION_FILE)
    @info "$SOLUTION_FILE not found, attempting to recreate one."

    SDP_problem, varP = PropertyT.SOS_problem_primal(elt, Δ, block_decomposition; upper_bound=upper_bound)

    with_SCS = JuMP.with_optimizer(SCS.Optimizer, linear_solver=SCS.DirectSolver,
                             max_iters=500_000,
                             eps=scseps,
                             alpha=1.5,
                             acceleration_lookback=0,
                             warm_start=true)

    let
        ws = nothing
        if isfile(WARMSTART_FILE)
            ws = load(WARMSTART_FILE, "warmstart")
        end

        status = JuMP.MOI.NO_SOLUTION
        while status != JuMP.MOI.OPTIMAL
            SOLVERLOG_FILE = joinpath(fullpath, "$(ELT_STRING)_solver_$(now()).log")
            @info "Recording solvers progress in" SOLVERLOG_FILE
            @time status, ws = PropertyT.solve(SOLVERLOG_FILE, SDP_problem, with_SCS, ws);
            λ = JuMP.value(SDP_problem[:λ])
            Ps = [JuMP.value.(P) for P in varP]

            if all((!isnan).(ws[1])) # solution looks valid
                @info "Warmstart looks valid, overriding $WARMSTART_FILE"
                save(WARMSTART_FILE,
                    "warmstart", (ws.primal, ws.dual, ws.slack), "Ps", Ps, "λ", λ)
                save(WARMSTART_FILE[1:end-4]*"$(now())"*".jld",
                    "warmstart", (ws.primal, ws.dual, ws.slack), "Ps", Ps, "λ", λ)
                begin
                    @info "Reconstructing Q..."
                    Qs = real.(sqrt.(Ps))
                    @time Q = PropertyT.reconstruct(Qs, block_decomposition);
                    save(SOLUTION_FILE, "λ", λ, "Q", Q)
                end
            else
                @warn "No valid solution was saved!" status
                break
            end
        end
    end
end

@info "Checking the sum of squares solution for Adj_N + $K·Op_N - $upper_bound·Δ_N"
Q, λ = load(SOLUTION_FILE, "Q", "λ")

function SOS_residual(eoi::GroupRingElem, Q::Matrix)
    RG = parent(eoi)
    @time sos = PropertyT.compute_SOS(RG, Q);
    return eoi - sos
end

let EOI = elt - λ*Δ
    residual = SOS_residual(EOI, Q)
    @info "In floating point arithmetic the ℓ₁-norm of the residual" norm(residual, 1);
end

let EOI_int = elt - @interval(λ)*Δ;
    Q_int, check_columns_in_augmentation_ideal = PropertyT.augIdproj(Interval, Q);
    @assert check_columns_in_augmentation_ideal

    residual_int = SOS_residual(EOI_int, Q_int);
    @info "In interval arithmetic the ℓ₁-norm of the residual" norm(residual_int, 1);

    λ_cert = @interval(λ) - 2*norm(residual_int,1)
    @info "λ is certified to be > " λ_cert.lo
    @info "i.e Adj_$N + $K·Op_$N - ($(λ_cert.lo))·Δ_$N ∈ Σ²₂ ISAut(F_$N)"
end


text:
"""
    GlyphCoordinates{T} <: AbstractVector{T}

A sparse representation of a [`Glyph`](@ref).
"""
struct GlyphCoordinates{T} <: AbstractVector{T}
    v::Vector{SVector{2,T}}
    sz::Tuple{Int,Int}
end

Base.size(g::GlyphCoordinates) = size(g.v)
Base.getindex(g::GlyphCoordinates, i::Int) = getindex(g.v, i)
Base.getindex(g::GlyphCoordinates, I...) = getindex(g.v, I...)
Base.IndexStyle(g::GlyphCoordinates) = IndexLinear()

GlyphCoordinates(args...) = GlyphCoordinates{Float64}(args...)

function GlyphCoordinates{T}(g::Glyph) where {T}
    GlyphCoordinates([SVector{2,T}(Tuple(ci))
                      for ci in CartesianIndices(g) if !iszero(g[ci])], size(g))
end

const COORDS_CACHE = Dict{Char,GlyphCoordinates{Float64}}()

function GlyphCoordinates{Float64}(c::Char)
    get!(COORDS_CACHE, c) do
        GlyphCoordinates{Float64}(Glyph(c))
    end
end

# fallback for generic types
GlyphCoordinates{T}(c::Char) where {T} = GlyphCoordinates{T}(Glyph(c))

# Use the character cache `COORDS_CACHE`
function GlyphCoordinates{T}(s::AbstractString) where {T}
    gcs = [GlyphCoordinates{T}(c) for c in s]
    L = sum(length, gcs)
    v = Vector{SVector{2,T}}(undef, L)
    shift = SVector{2,Int}(0, 0)
    j = 1
    for gc in gcs
        l = length(gc.v)
        v[j:j+l-1] .= gc.v .+ Ref(shift)
        j += l
        shift += SVector{2,Int}(0, gc.sz[2])
    end
    GlyphCoordinates{T}(v, Tuple(shift + SVector(16, 0)))
end

function printglyph(io, g::GlyphCoordinates{T}; symbols=("#", " ")) where {T}
    for r = 1:g.sz[1]
        for c = 1:g.sz[2]
            if SVector{2,T}(r, c) ∈ g.v
                print(io, symbols[1])
            else
                print(io, symbols[2])
            end
        end
        println(io)
    end
end


text:
function run_import_pkgs(argv=ARGS)

    argset = ArgParse.ArgParseSettings()
    ArgParse.@add_arg_table! argset begin
        "rootdir"
            help="Root directory"
            arg_type = String
            default = pwd()
    end

    parsed_args = ArgParse.parse_args(argv, argset)
    rootdir = parsed_args["rootdir"]

    _walk_pkgs(;rootdir) do path, proj

        ispackage = haskey(proj, "name")
        if ispackage
            name = proj["name"]
            try;
                cmd = _Cmd(["julia", "--startup-file=no", "--project=$(path)", "-e",
                    "import $(name); println(\"$(name)\", \" imported\")"])
                run(cmd; wait = true)
            catch err;
                @error("Failing importing ", name, err)
            end
        else
        end
    end
end


text:
module InverseGaussTransform

using StokesParticles
using Plots

include("ctd_functions.jl")

end


text:
using DataFrames, QuasiCopula, LinearAlgebra, GLM, RDatasets
using Test

Mmmec = dataset("mlmRev", "Mmmec");
df = Mmmec
y = :Deaths
grouping = :Region
covariates = [:UVB]

### Now with two variance components V1 = ones(di, di), V2 = Identity
function make_Vs(
    df::DataFrame,
    y::Symbol,
    grouping::Symbol)
    groups = unique(df[!, grouping])
    n = length(groups)
    V = Vector{Vector{Matrix{Float64}}}(undef, n)
    for (i, grp) in enumerate(groups)
        gidx = df[!, grouping] .== grp
        ni = count(gidx)
        V[i] = [ones(ni, ni), Matrix(I, ni, ni)]
    end
    V
end

V = make_Vs(df, y, grouping);

#### POISSON ####
d = Poisson()
link = LogLink()
Poisson_VC_model = VC_model(df, y, grouping, covariates, d, link)

@test typeof(Poisson_VC_model) == GLMCopulaVCModel{Float64, Poisson{Float64}, LogLink}
# test that we have all ones for random intercept covariance structure by default
@test Poisson_VC_model.data[1].V[1][1] == 1.0

# two vc
Poisson_2VC_model = VC_model(df, y, grouping, covariates, V, d, link)
@test typeof(Poisson_2VC_model) == GLMCopulaVCModel{Float64, Poisson{Float64}, LogLink}
@test length(Poisson_2VC_model.data[1].V) == 2

# intercept only model
Poisson_VC_model_intercept_only = VC_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(Poisson_VC_model_intercept_only.data[1].X) == Poisson_VC_model_intercept_only.data[1].n

# intercept only model with V and penalty
Poisson_2VC_model_intercept_only = VC_model(df, y, grouping, V, d, link; penalized = true)
# check if default to intercept only
@test Poisson_2VC_model_intercept_only.penalized == true
@test sum(Poisson_2VC_model_intercept_only.data[1].X) == Poisson_2VC_model_intercept_only.data[1].n

#### nb ####
d = NegativeBinomial()
link = LogLink()
nb_VC_model = VC_model(df, y, grouping, covariates, d, link)

@test typeof(nb_VC_model) == NBCopulaVCModel{Float64, NegativeBinomial{Float64}, LogLink}
# test that we have all ones for random intercept covariance structure by default
@test nb_VC_model.data[1].V[1][1] == 1.0

# two vc
nb_2VC_model = VC_model(df, y, grouping, covariates, V, d, link)
@test typeof(nb_2VC_model) == NBCopulaVCModel{Float64, NegativeBinomial{Float64}, LogLink}
@test length(nb_2VC_model.data[1].V) == 2

# intercept only model
nb_VC_model_intercept_only = VC_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(nb_VC_model_intercept_only.data[1].X) == nb_VC_model_intercept_only.data[1].n

# intercept only model with V and penalty
nb_2VC_model_intercept_only = VC_model(df, y, grouping, V, d, link; penalized = true)
# check if default to intercept only
@test nb_2VC_model_intercept_only.penalized == true
@test sum(nb_2VC_model_intercept_only.data[1].X) == nb_2VC_model_intercept_only.data[1].n

#### gaussian ####
d = Normal()
link = IdentityLink()
gaussian_VC_model = VC_model(df, y, grouping, covariates, d, link)

@test typeof(gaussian_VC_model) == GaussianCopulaVCModel{Float64}
# test that we have all ones for random intercept covariance structure by default
@test gaussian_VC_model.data[1].V[1][1] == 1.0

# two vc
gaussian_2VC_model = VC_model(df, y, grouping, covariates, V, d, link)
@test typeof(gaussian_2VC_model) == GaussianCopulaVCModel{Float64}
@test length(gaussian_2VC_model.data[1].V) == 2

# intercept only model
gaussian_VC_model_intercept_only = VC_model(df, y, grouping, d, link)
# check if default to intercept only
@test sum(gaussian_VC_model_intercept_only.data[1].X) == gaussian_VC_model_intercept_only.data[1].n

# intercept only model with V and penalty
gaussian_2VC_model_intercept_only = VC_model(df, y, grouping, V, d, link; penalized = true)
# check if default to intercept only
@test gaussian_2VC_model_intercept_only.penalized == true
@test sum(gaussian_2VC_model_intercept_only.data[1].X) == gaussian_2VC_model_intercept_only.data[1].n


text:
input = 347991

# Part 1
function side_length(number:: Int)
  root = ceil(sqrt(number))
  if root % 2 == 0
    return root + 1
  else
    return root
  end
end

function get_position(number:: Int)
  slength = side_length(number)
  x = floor(slength / 2)

  side_centers = slength^2 - ((slength - 1) * collect(1:4)) - floor(slength / 2)
  y = minimum(abs.(side_centers - number))

  return x, y
end

(x, y) = abs.(get_position(input))
println("Part 1: $(x + y) steps from ($x, $y)")

# Part 2
mutable struct NegativeIndexSquareArray
  size::Int
  content::Array{Int, 3}

  NegativeIndexSquareArray(size) = new(size, zeros(size, size, 4))
end

function Base.getindex(array::NegativeIndexSquareArray, x, y)
  sx = sign(x) == 0? 1 : sign(x)
  sy = sign(y) == 0? 1 : sign(y)
  (x, y) = abs.([x, y])

  if sx == 1 && sy == 1
    array.content[x + 1, y + 1, 1]
  elseif sx == 1 && sy == -1
    array.content[x + 1, y, 2]
  elseif sx == -1 && sy == 1
    array.content[x, y + 1, 3]
  else # Both negative
    array.content[x, y, 4]
  end
end

function Base.setindex!(array::NegativeIndexSquareArray, value, x, y)
  sx = sign(x) == 0? 1 : sign(x)
  sy = sign(y) == 0? 1 : sign(y)
  (x, y) = abs.([x, y])

  if sx == 1 && sy == 1
    array.content[x + 1, y + 1, 1] = value
  elseif sx == 1 && sy == -1
    array.content[x + 1, y, 2] = value
  elseif sx == -1 && sy == 1
    array.content[x, y + 1, 3] = value
  else # Both negative
    array.content[x, y, 4] = value
  end
end

spiral = NegativeIndexSquareArray(10)
spiral[0, 0] = 1
(x, y) = [1, 0]
direction = im

while true
  spiral[x, y] = spiral[x - 1, y] + spiral[x - 1, y - 1] +
    spiral[x, y - 1] + spiral[x + 1, y - 1] + spiral[x + 1, y] +
    spiral[x + 1, y + 1] + spiral[x, y + 1] + spiral[x - 1, y + 1]

  if spiral[x, y] > input
    println("Part 2: $(spiral[x, y]) at ($x, $y)")
    break
  end

  left = direction * im
  if spiral[x + real(left), y + imag(left)] == 0
    direction = left
  end
  x += real(direction)
  y += imag(direction)
end


text:
using Test, LinearAlgebra, InvertibleNetworks, Flux

# Input
nx = 64
ny = 64
n_in = 10
n_hidden = 20
batchsize = 2

X = glorot_uniform(nx, ny, n_in, batchsize)

# Flux network
model = Chain(
    Conv((3,3), n_in => n_hidden; pad=1),
    BatchNorm(n_hidden, relu),
    Conv((3,3), n_hidden => n_hidden; pad=1),
    BatchNorm(n_hidden, relu),
    Conv((3,3), n_hidden => n_in; pad=1),
    BatchNorm(n_in, relu)
)

# Create Flux CNN block
FB = FluxBlock(model)

# Evaluate forward
Y = FB.forward(X)

# Backpropagate data residual
ΔY = Y.*.1f0  # data residual
ΔX = FB.backward(ΔY, X)

# Update weights
Params = get_params(FB)
α = 1f-2
for p in Params
    p.data[:] += α*p.grad[:]    # need to set values via [:]
end
clear_grad!(FB)


text:
module TestWiringDiagramAlgebras
using Test

using Tables, TypedTables

using Catlab.CategoricalAlgebra, Catlab.CategoricalAlgebra.FinSets
using Catlab.Graphs, Catlab.WiringDiagrams, Catlab.Programs.RelationalPrograms

tuples(args...) = sort!(collect(zip(args...)))

# UWD algebra of structured multicospans
########################################

const OpenGraphOb, OpenGraph = OpenCSetTypes(Graph, :V)

# Directed fork as open graph with one input and two outputs.
g0 = Graph(4)
add_edges!(g0, [1,2,2], [2,3,4])
g = OpenGraph(g0, FinFunction([1],4), FinFunction([3,4],4))

# Opposite of previous graph.
h0 = Graph(4)
add_edges!(h0, [1,2,3], [3,3,4])
h = OpenGraph(h0, FinFunction([1,2],4), FinFunction([4],4))

# Sequential composition.
seq = @relation (a,c) where (a,b,c) begin
  g(a,b)
  h(b,c)
end
k = oapply(seq, Dict(:g => g, :h => h))
@test length(legs(k)) == 2
@test feet(k) == [first(feet(g)), last(feet(h))]
k0 = apex(k)
@test (nv(k0), ne(k0)) == (6, 6)
@test tuples(src(k0), tgt(k0)) == [(1,2), (2,3), (2,4), (3,5), (4,5), (5,6)]
@test [ collect(leg[:V]) for leg in legs(k) ] == [[1], [6]]

# Parallel composition.
para = @relation (a,c,b,d) where (a,b,c,d) begin
  g(a,b)
  h(c,d)
end
k = oapply(para, Dict(:g => g, :h => h))
@test length(legs(k)) == 4
@test feet(k) == [first(feet(g)), first(feet(h)), last(feet(g)), last(feet(h))]
k0 = apex(k)
@test (nv(k0), ne(k0)) == (8, 6)
@test tuples(src(k0), tgt(k0)) == [(1,2), (2,3), (2,4), (5,7), (6,7), (7,8)]
@test [ collect(leg[:V]) for leg in legs(k) ] == [[1], [5,6], [3,4], [8]]

# Identity for sequential composition.
seq_id = @relation (a,a) where (a,) begin end
k = oapply(seq_id, Dict{Symbol,OpenGraph}(), Dict(:a => FinSet(3)))
@test apex(k) == Graph(3)
@test feet(k) == [FinSet(3), FinSet(3)]

# Queries of ACSets
###################

paths2 = @relation (start=start, stop=stop) begin
  E(src=start, tgt=mid)
  E(src=mid, tgt=stop)
end

# Graph underlying a commutative squares.
square = Graph(4)
add_edges!(square, [1,1,2,3], [2,3,4,4])
result = query(square, paths2)
@test result == Table((start=[1,1], stop=[4,4]))

# Graph underlying a pasting of two commutative squares.
squares2 = copy(square)
add_vertices!(squares2, 2)
add_edges!(squares2, [2,4,5], [5,6,6])
result = query(squares2, paths2)
@test tuples(columns(result)...) == [(1,4), (1,4), (1,5), (2,6), (2,6), (3,6)]

result = query(squares2, paths2, (start=1,))
@test tuples(columns(result)...) == [(1,4), (1,4), (1,5)]
result = query(squares2, paths2, (start=1, stop=4))
@test result == Table((start=[1,1], stop=[4,4]))

cycles3 = @relation (edge1=e, edge2=f, edge3=g) where (e,f,g,u,v,w) begin
  E(_id=e, src=u, tgt=v)
  E(_id=f, src=v, tgt=w)
  E(_id=g, src=w, tgt=u)
end

g = cycle_graph(Graph, 3)
result = query(g, cycles3)
@test tuples(columns(result)...) == [(1,2,3), (2,3,1), (3,1,2)]
result = query(g, cycles3, (v=1,))
@test result == Table((edge1=[3], edge2=[1], edge3=[2]))

@test isempty(query(cycle_graph(Graph, 4), cycles3))

end


text:
module FFD

export AbstractMappingType, Mapping, PumiMapping
export PumiBoundingBox, calcKnot, controlPoint, calcParametricMappingLinear
export calcParametricMappingNonlinear, evalVolume, evalSurface
export writeControlPointsVTS, evaldXdControlPointProduct, evaldXdControlPointTransposeProduct

export initializeFFD, commitToPumi

export numLinearPlaneConstraints, countVarsLinearPlaneConstraints!, setLinearPlaneConstraints!
export numLinearCornerConstraints, countVarsLinearCornerConstraints!, setLinearCornerConstraints!
export numLinearStretchConstraints, countVarsLinearStretchConstraints!, setLinearStretchConstraints!
export numLinearRootConstraints, countVarsLinearRootConstraints!, setLinearRootConstraints!

export getControlPoints, setControlPoints

#push!(LOAD_PATH, joinpath(Pkg.dir("PumiInterface"), "src"))

using ArrayViews
using PumiInterface
using PdePumiInterface
using ODLCommonTools
using SummationByParts
using WriteVTK

using MPI

# initialize MPI, and arrange for its finalization if so
function finalizeMPI()
#  println("running atexit hook for MPI")
  if MPI.Initialized()
#    println("finalizing MPI")
    MPI.Finalize()
  end
end

if !MPI.Initialized()
  MPI.Init()
  atexit(finalizeMPI)
end

# Abstract Type definition
abstract type AbstractMappingType{Tffd} <: Any end # Abstract Mapping type for creating a different Mapping type
abstract type AbstractBoundingBox{Tffd} <: Any end

@doc """
### Mapping

A type for creating mapping objects. The mapping uses a uniform knot
distribution along the 3 dimensions in the parametric space.  An object of this
type can be defined by calling an inner constructor with the following arguments
in sequence

* Number of dimensions (2 or 3)
* Array of order of B-splines in all dimensions
* Number of control points along every direction
* Number of nodes (embedded geometry points) along every direction

**Method 2 constructor**
  function PumiMapping(ndim::Int, order::AbstractArray{Int,1},
                       nctl::AbstractArray{Int,1}, mesh_info::AbstractArray{Int,1};
                       full_geom=true, bc_nums=[0])

*  `ndim` : Number of dimensions
*  `order`: Aray of B-spline order in the 3 dimenstions
*  `nctl` : number of control points along the 3 dimenstions
*  `mesh_info` : information about the geometry being embedded. If the entire
                 geometry is being embedded then
```
  mesh_info = [num_nodes_per_element, num_elements]
```
                 else,
```
  mesh_info = [sbp.facenodes, num_boundary_faces, num_geometric_faces]
```
*  `full_geom` : Bool. `true` for embedding the entire geometry. `false` if you
                 want to embedded only certain faces of a geometry
*  `bc_nums` : Array of boundary condition numbers that identify the surface
               being embedded in the FFD box

**Members**

*  `ndim`    : Number of dimensions (2D or 3D)
*  `nctl`    : Array of number of control points in each direction
*  `numnodes`: Array of number of nodes in each direction
*  `order`   : Array of order of B-splines in each direction
*  `xi`      : Array of parametric coordinates of the nodes of the geometry
               embedded within. The parametric axes within FFD are called
               (s,t,u)
*  `cp_xyz`  : Control point coordinates in the (x,y,z) space
*  `edge_knot` : 3D Array of edge knot vectors. dim1 = Knot vector, dim2 = edge
                 number (1:4), dim3 = direction (di)
*  `edge_param`: 3D array to store edge spacing parameters. dim1 = A or b,
                 dim2 = edge number (1:4), dim3 = direction (di)
*  `aj`   :
*  `dl` and `dr`   : Working array for computing spline basis and its derivatives
*  `work` :

"""->

mutable struct Mapping{Tffd} <: AbstractMappingType{Tffd}

  ndim::Int                     # Mapping object to indicate 2D or 3D
  nctl::AbstractArray{Int, 1}   # Number of control points in each of the 3 dimensions
  numnodes::AbstractArray{Int, 1} # Number of nodes in each direction
  order::AbstractArray{Int, 1}  # Order of B-spline in each direction

  xi::AbstractArray{AbstractFloat, 4}     # Coordinate values
  cp_xyz::AbstractArray{AbstractFloat, 4} # Cartesian coordinates of control points
  edge_knot::AbstractArray{Vector{AbstractFloat}, 1}  # edge knot vectors

  # Working arrays
  aj::AbstractArray{AbstractFloat, 3}
  dl::AbstractArray{AbstractFloat, 2}
  dr::AbstractArray{AbstractFloat, 2}
  work::AbstractArray{AbstractFloat, 4}

  evalVolume::Function

  function Mapping{Tffd}(dim, k, ncpts, nnodes) where {Tffd}

    # Assertion statements to prevent errors
    @assert dim >= 2 "Only 2D and 3D valid"
    for i = 1:3
      @assert ncpts[i] > 0 "Number of control points specified is <= 0 in $i direction"
      @assert nnodes[i] > 0 "Number of nodes specified is <= 0 in $i direction"
      @assert k[i] > 0 "Order of B-spline specified is <= 0 in $i direction"
    end

    # Define max_wrk = number of work elements at each control point
    const max_work = 2*6  # 2*n_variables

    ndim = dim  # To indicate a 3D Mapping object is being created
    nctl = zeros(Int, 3)
    numnodes = zeros(Int, 3)
    order = zeros(Int, 3)

    nctl[:] = ncpts[:]    # Set number of control points
    numnodes[:] = nnodes[:] # Set map refinement level in each coordinate direction
    order[:] = k[:]       # Set the order of B-splines

    for i = 1:3
      @assert nctl[i] > 0
      @assert numnodes[i] > 0
      @assert order[i] > 0
    end

    # Allocate and initialize mapping arrays
    max_order = maximum(order)  # Highest order among 3 dimensions
    max_knot = max_order + maximum(nctl) # Maximum number of knots among 3 dimensions

    xi = zeros(numnodes[1], numnodes[2], numnodes[3], 3)
    cp_xyz = zeros(3,nctl[1], nctl[2], nctl[3])
    edge_knot = Array{Vector{AbstractFloat}}(3)
    for i = 1:3
      edge_knot[i] = zeros(AbstractFloat, nctl[i]+order[i])
    end

    aj = zeros(3, max_order, 3)
    dl = zeros(max_order-1, 3)
    dr = zeros(max_order-1, 3)
    work = zeros(nctl[1], nctl[2], nctl[3], max_work)

    new(ndim, nctl, numnodes, order, xi, cp_xyz, edge_knot, aj, dl, dr, work)

  end  # End constructor

end  # End Mapping

@doc """
### PumiMapping

  This subtype of [`AbstractMapping`](@ref) enables the use of Free Form
  Deformation with a Pumi mesh.

  Users should *not* call this function directly, they should call
  [`initializeFFD`](@ref) instead.

  **Public Fields**

   * ndim: dimensionality of the mesh
   * cp_xyz: a 3 x `nctl[1]` x `nctl[2]` x `nctl[3]` array containing the
             x,y, and z coordinates of each control point.  Note that
             the control point grid is always 3 dimensional, even when the
             mesh is 2 dimensional.  This array is read-only
    * n_face: a Pumi apf::Numbering object that numbers the face nodes
    * numFacePts: the number of points on the face
    * nctl: array of length 3, containing the number of control points in each
            direction.  For 2D meshes, the number of control points in the
            z direction must be 2.

  **Private Fields*

   * bc_nums: the array of boundary condition numbers that define the
              surface enveloped by the FFD control points
   * face_verts: vector of apf::MeshEntity* (Ptr{Void}) containing the
                 vertices in the order defined by `n_face`.  Length
                 `numFacePts`
   * map_cs: a PumiMapping{Complex128}, used for doing complex step when
             Tffd is real

  **Constructors**

   ***Constructor 1***

    This zero argument constructor produces an object with all fields
    zero.  Useful for constructing a dummy object

   ***Constructor 2***

    This constructor is the main constructor.  Its argument are:

     * ndims: number of dimensions (must be the same as mesh.dim
     * order: array of length 3 specifying the degree of the B-splines in
              each direction
     * nctl: array of length 3 specifying the number of control points in
             each direction
     * mesh: the Pumi mesh object
     * n_face: a Pumi apf::Numbering* hat numbers the surface points, default
               to C_NULL.  If unspecified, a new surface numbering will be
               created.
     * numFacePts: number of face points numbered by n_face.  Defaults to
                   zero.  This argument should only be specified in `n_face` is
                   as well.
     * bc_nums: array of integers specifying the boundary condition numbers
                (used in the construction of the `mesh`) that identify the
                surface to envelop in the FFD volume.  This is a keyword
                argument

  ***Constructor 3***

  Copy constructor.  Allows copying a PumiMapping.  The new one can have a
  different type parameter `Tffd`.  Arguments:

   * map: a PumiMapping
"""->

mutable struct PumiMapping{Tffd} <: AbstractMappingType{Tffd}

  ndim::Int                     # Mapping object to indicate 2D or 3D
  nctl::Array{Int, 1}   # Number of control points in each of the 3 dimensions
  order::Array{Int, 1}  # Order of B-spline in each direction

  xi::Array       # Paramaetric coordinates of input geometry
  cp_xyz::ROView{Tffd, 4, Array{Tffd, 4}}
  _cp_xyz::Array{Tffd, 4} # Cartesian coordinates of control points
  edge_knot::Array{Vector{Tffd}, 1}  # edge knot vectors
  bc_nums::Array{Int,1}
  cp_idx::Array{Int,4}  # index assigned to each CP coordinate

  # Working arrays
  aj::Array{Tffd, 3}
  dl::Array{Tffd, 2}
  dr::Array{Tffd, 2}
  work::Array{Tffd, 4}

  n_face::Ptr{Void}  # apf::Numbering*, numbering the face points
  numFacePts::Int  # number of face points
  face_verts::Array{Ptr{Void}, 1}  # face points, in order

  map_cs::PumiMapping{Complex128}  # map for complex stepping even when Tffd
                                   # is real
  vertices::Array{Tffd, 2}  # dim x numFacePts array, used by map_cs

  function PumiMapping{Tffd}(ndim::Int, order::AbstractArray{Int,1},
                       nctl::AbstractArray{Int,1}, mesh::AbstractMesh,
                       n_face::Ptr{Void}=C_NULL, numFacePts=0;
                       bc_nums::AbstractArray{Int,1}=[0]) where {Tffd}

    #TODO: remove dims, use mesh.dim instead
    map = new()
    # Check if the input arguments are valid
    @assert ndim >= 2 "Only 2D and 3D valid"
    @assert mesh.coord_order <= 2
    for i = 1:3
      @assert order[i] > 0 "Order cannot be 0 in the $i direction"
      @assert nctl[i] > 0 "number of control points cannot be 0 in $i direction"
    end

    # Define max_wrk = number of work elements at each control point
    const max_work = 2*6  # 2*n_variables

    map.ndim = mesh.dim
    map.order = order
    map.nctl = nctl
    map._cp_xyz = zeros(3, nctl[1], nctl[2], nctl[3])
    map.cp_xyz = ROView(map._cp_xyz)

    map.edge_knot = Array{Vector{Tffd}}(3)
    for i = 1:3
      map.edge_knot[i] = zeros(Tffd, nctl[i]+order[i])
    end

    # if the user didn't provide a surface point numbering, make one ourselves
    if n_face == C_NULL
      # call PdePumiInterface function
      numFacePts, n_face, face_verts = numberSurfacePoints(mesh, bc_nums)
    else  # get the face_verts array
      face_verts = Array{Ptr{Void}}(numFacePts)
      for vert in mesh.verts
        n_v = getNumberJ(n_face, vert, 0, 0)
        if n_v <= numFacePts
          face_verts[n_v] = vert
        end
      end

      # get edges
      # doing all the edges after the vertices results in a non-cache
      # optimal order
      if hasNodesIn(mesh.coordshape_ptr, 1)
        for edge in mesh.edges
          n_e = getNumberJ(n_face, edge, 0, 0)
          if n_e <= numFacePts
            face_verts[n_e] = edge
          end
        end
      end
    end  # end if n_face = C_NULL
    map.n_face = n_face
    map.numFacePts = numFacePts
    map.face_verts = face_verts
    map.xi = zeros(Tffd, 3, numFacePts)
    map.bc_nums = bc_nums

    # Allocate and initialize mapping arrays
    max_order = maximum(order)  # Highest order among 3 dimensions
    max_knot = max_order + maximum(nctl) # Maximum number of knots among 3 dimensions

    # use simple logical indexing for CP coordinates
    map.cp_idx = zeros(3, nctl[1], nctl[2], nctl[3])
    ptr = 1
    for k = 1:nctl[3]
      for j = 1:nctl[2]
        for i = 1:nctl[1]
          for di = 1:3
            map.cp_idx[di,i,j,k] = ptr
            ptr += 1
          end
        end
      end
    end

    map.aj = zeros(3, max_order, 3)
    map.dl = zeros(max_order-1, 3)
    map.dr = zeros(max_order-1, 3)
    map.work = zeros(max_work, nctl[1], nctl[2], nctl[3])

    # do this in initializeFFD because map isn't fully initialized yet
#    map.map_cs = PumiMapping{Complex128}(map)
    map.vertices = Array{Tffd}(0, 0)

    return map
  end

  function PumiMapping{Tffd}()  where {Tffd}
    map = new()
    map.ndim = 0
    map.nctl = Int[]
    map.order = Int[]
    map.xi = zeros(Tffd, 0, 0)
    map._cp_xyz = zeros(Tffd, 0, 0, 0, 0)
    map.cp_xyz = ROView(map._cp_xyz)
    map.edge_knot = Array{Vector{Tffd}}(0)
    map.bc_nums = Int[]
    map.cp_idx = zeros(Tffd, 0, 0, 0, 0)
    map.aj = zeros(Tffd, 0, 0, 0)
    map.dl = zeros(Tffd, 0, 0)
    map.dr = zeros(Tffd, 0, 0)
    map.work = zeros(Tffd, 0, 0, 0, 0)

    map.n_face = C_NULL
    map.numFacePts = 0
    map.face_verts = Ptr{Void}[]
    map.map_cs = PumiMapping{Tffd}(map)
    map.vertices = Array{Tffd}(0, 0)

    return map
  end

  # constructs a new mapping from a fully initialized one, with possibly
  # different parameter Tffd
  function PumiMapping{Tffd}(map_old::PumiMapping)  where {Tffd}

    map = new()
    map.ndim = map_old.ndim
    map.nctl = copy(map_old.nctl)
    map.order = copy(map_old.order)
    map.xi = zeros(Tffd, size(map_old.xi)); copy!(map.xi, map_old.xi)
    map._cp_xyz = zeros(Tffd, size(map_old.cp_xyz)); copy!(map._cp_xyz, map_old._cp_xyz)
    map.cp_xyz = ROView(map._cp_xyz)
    map.edge_knot = Array{Vector{Tffd}}(length(map_old.edge_knot))
    for i=1:length(map.edge_knot)
      map.edge_knot[i] = copy(map_old.edge_knot[i])
    end
    map.bc_nums = copy(map_old.bc_nums)
    map.cp_idx = copy(map_old.cp_idx)

    map.aj = zeros(Tffd, size(map_old.aj)); copy!(map.aj, map_old.aj)
    map.dl = zeros(Tffd, size(map_old.dl)); copy!(map.dl, map_old.dl)
    map.dr = zeros(Tffd, size(map_old.dr)); copy!(map.dr, map_old.dr)
    map.work = zeros(Tffd, size(map_old.work)); copy!(map.work, map_old.work)

    map.n_face = map_old.n_face
    map.numFacePts = map_old.numFacePts
    map.face_verts = map_old.face_verts  # reference, don't copy for this one
                                         # because it is large
    # leave map_cs uninitialized

    map.vertices = Array{Tffd}(map_old.ndim, map.numFacePts)
    return map
  end

end  # End type PumiMapping

include("knot.jl")
include("bounding_box.jl")
include("mapping_functions.jl")
include("control_point.jl")
include("span.jl")
include("b-splines.jl")
include("evaluations.jl")
include("constraints.jl")
include("pumi_specific_functions.jl")

@doc """
Routine to be called externally for initializing FreeFormDeformation

**Input**

* `mesh` : Pumi mesh
* `sbp`  : Summation-By-Parts operator
* `order`: Array of B-splone order along the 3 parametric directions (ξ, η, ζ)
* `nControlPts` : Array of nu,ber of control points along 3 parametric
                  directions (ξ, η, ζ)
* `offset` : FFD Bounding box offset form the embedded geometry along the 3
             physical directions (x, y, z)
* `bc_nums`: (Optional Argument) If partial geometry embedded, give the
                face/edge numbers of the embedded gemetry
* `n_face`: an apf::Numbering* (Ptr{Void}) numbering the face nodes, if
            not supplied a new numbering will be created
* `numFacePts`: the number of face points on the surface defined by
                `bc_nums` and numbered by `n_face`

"""->

function initializeFFD(mesh::AbstractMesh{Tmsh}, sbp::AbstractSBP,
                 order::AbstractArray{Int,1},
                 nControlPts::AbstractArray{Int,1},
                 offset::AbstractArray{Float64,1},
                 bc_nums::AbstractArray{Int,1}=[0],
                 n_face::Ptr{Void}=C_NULL,
                 numFacePts::Integer=0) where Tmsh

  # Create Mapping object
  ndim = mesh.dim
  map = PumiMapping{Tmsh}(ndim, order, nControlPts, mesh, n_face, numFacePts,
                              bc_nums=bc_nums)

  # Create knot vector
  calcKnot(map)

  # Create Bounding box
  # this works for curvilinear
  ffd_box = PumiBoundingBox{Tmsh}(map, mesh, sbp, offset)

  # Control points
  controlPoint(map, ffd_box)

  # Populate map.xi
  calcParametricMappingNonlinear(map, ffd_box, mesh, bc_nums)

  # now that map is fully initialized, copy it for the complex step calculations
  map.map_cs = PumiMapping{Complex128}(map)
  return map, ffd_box
end

"""
  This function copies the supplied array into map._cp_xyz, without checking
  the dimensionality first.

  **Inputs**

   * map: a PumiMapping object
   * cp_xyz: array, same size as map.cp_xyz, with new control point coordinates.
"""
function _setControlPoints(map::PumiMapping, cp_xyz::AbstractArray{T, 4}) where T

  for i=1:4
    @assert size(cp_xyz, i) == size(map._cp_xyz, i)
  end

  for i=1:length(cp_xyz)
    map._cp_xyz[i] = cp_xyz[i]
  end

  return nothing
end

"""
  This function updates the control point locations.  Users should not modify
  map.cp_xyz directly, they should use this function instead.  Derivative
  values will be incorrect if map.cp_xyz is modified.

  **Inputs**

   * map: a PumiMapping
   * cp_xyz: array with new control point locations.  In 2D, this array
             should be 2 x map.nctl[1] x map.nctl[2].  In 3D it should be
             3 x map.nctl[1] x map.nctl[2] x map.nctl[3]

  **Implementaton Notes**

  The internals of FFD always assume a 3 dimensional mesh, but for running
  2D problems having 2 control points along the z axis doesn't make sense.
  In 2D, this function updates both planes of control points simultaneously.
"""
function setControlPoints(map::PumiMapping, cp_xyz::AbstractArray{T, 4}) where T

  @assert map.ndim == 3

  # copy everything over
  _setControlPoints(map, cp_xyz)

  return nothing
end

function setControlPoints(map::PumiMapping, cp_xyz::AbstractArray{T, 3}) where T

  @assert map.ndim == 2
  @assert size(cp_xyz, 1) == 2
  @assert size(cp_xyz, 2) == size(map.cp_xyz, 2)
  @assert size(cp_xyz, 3) == size(map.cp_xyz, 3)
  @assert size(map.cp_xyz, 4) == 2

  # set cp_xyz for both z = zmin and z = zmax
  # note that the derivative calculations need to take this into account

  for i=1:size(cp_xyz, 3)
    for j=1:size(cp_xyz, 2)
      for k=1:size(cp_xyz, 1)
        map._cp_xyz[k, j, i, 1] = cp_xyz[k, j, i]
        map._cp_xyz[k, j, i, 2] = cp_xyz[k, j, i]
      end
    end
  end

  return nothing
end

"""
  This function returns a copy of the current control point locations, suitable
  to be passed into [`setControlPoints`](@ref)  See that function for details

  Note that this function is type unstable.  Users are generally better off
  allocating their own array.

  **Inputs**

   * map: PumiMapping

  **Outputs**

   * array of control point locations
"""
function getControlPoints(map::PumiMapping)

  if map.ndim == 2
    return map._cp_xyz[1:2, :, :, 1]
  else
    return copy(map._cp_xyz)
  end

  return nothing
end

"""
  Second method for getting control points, where the user supplied the array.
  The requires size of the array is described by [`setControlPoints`](@ref)

  **Inputs**

   * map: PumiMapping

  **Input/Outputs**

   * cp_xyz: array to be overwritten with control point coordinates
"""
function getControlPoints(map::PumiMapping, cp_xyz::AbstractArray{T, 3}) where T

  @assert map.ndim == 2
  @assert size(cp_xyz, 1) == 2
  @assert size(cp_xyz, 2) == map.nctl[1]
  @assert size(cp_xyz, 3) == map.nctl[2]

  for i=1:size(cp_xyz, 3)
    for j=1:size(cp_xyz, 2)
      for k=1:size(cp_xyz, 1)
        cp_xyz[k, j, i] = map._cp_xyz[k, j, i, 1]
      end
    end
  end

  return nothing
end

function getControlPoints(map::PumiMapping, cp_xyz::AbstractArray{T, 4}) where T

  @assert map.ndim == 3
  @assert size(cp_xyz, 1) == 3
  @assert size(cp_xyz, 2) == map.nctl[1]
  @assert size(cp_xyz, 3) == map.nctl[2]
  @assert size(cp_xyz, 3) == map.nctl[3]

  copy!(cp_xyz, map._cp_xyz)

  return nothing
end

"""
  This function checks if an array formatted similarly to map.cp_xyz
  has the same coordinates for the z = zmin and z = zmax.  This is the
  condition enforced by [`setControlPoints`](@ref) for 2D arrays.

  **Inputs**

   * cp_xyz: array like cp_xyz
"""
function check2DSymmetry(cp_xyz::AbstractArray{T, 4}, tol=1e-13) where T

  @assert size(cp_xyz, 4) == 2

  nonsym = false
  for i=1:size(cp_xyz, 3)
    for j=1:size(cp_xyz, 2)
      for k=1:2
        flag = abs(cp_xyz[k, j, i, 1] - cp_xyz[k, j, i, 2]) > tol
        nonsym = nonsym || flag
      end
    end
  end

  @assert !nonsym "cp_xyz is non-symmetric about the x-y plane.  Did you modify PumiMapping.cp_xyz directly?"

  return nothing
end

#=
@doc """
### FreeFormDeformation.defineMapXi

"""->

function defineMapXi(mesh::AbstractMesh, bc_nums::AbstractArray{Int,1},
                     xi::AbstractArray)

  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums =sview(mesh.bndryfaces, idx_range)
    nfaces = length(bndry_facenums)
    xi[idx] = zeros(3,mesh.coord_numNodesPerFace,nfaces)
  end

  return nothing
end
=#

@doc """
###FreeFormDeformation.defineVertices

This function defines the shape of the vertices array that is used to update
the Pumi mesh after FFD when a geometric face is paramtereized.

**Arguments**

* `mesh` : Pumi DG mesh
* `bc_nums` : Array of boundary condition numbers over which the number
                 of element faces needs to be computed
* `vertices` : Array of arrays holding the coodinates of the updated vertices
               shape = vertices[n_bc_nums][mesh.dim, mesh.coord_numNodesPerFace, n_elem_faces]

"""->

function defineVertices(mesh::AbstractDGMesh{Tmsh}, bc_nums::AbstractArray{Int,1},
                  vertices::AbstractArray) where Tmsh

  for (idx, itr) in enumerate(bc_nums)
    start_index = mesh.bndry_offsets[itr]
    end_index = mesh.bndry_offsets[itr+1]
    idx_range = start_index:(end_index-1)
    bndry_facenums = sview(mesh.bndryfaces, idx_range)
    nfaces = length(bndry_facenums)
    arr_dim = mesh.dim
    vertices[idx] = zeros(Tmsh, arr_dim, mesh.coord_numNodesPerFace, nfaces)
    for i = 1:nfaces
      bndry_i = bndry_facenums[i]
      for j=1:mesh.coord_numNodesPerFace
        v_j = mesh.topo.face_verts[j, bndry_i.face]
        vertices[idx][1:mesh.dim,j,i] = mesh.vert_coords[1:mesh.dim,v_j,bndry_i.element]
      end
    end # End for i = 1:nfaces
  end

  return nothing
end

#=
@doc """
###FreeFormDeformation.commitToPumi

Uses the modified vertex coordinates to update the Pumi mesh

"""->

function commitToPumi{Tffd, Tmsh}(map::PumiMapping{Tffd},
                      mesh::AbstractDGMesh{Tmsh}, sbp, vertices, opts)

  if map.full_geom == true
    println("Tmsh = $Tmsh")
    for i = 1:mesh.numEl
      update_coords(mesh, i, real(vertices[:,:,i]))
    end
  else
    for itr = 1:length(map.geom_faces)
      geom_face_number = map.geom_faces[itr]
      # get the boundary array associated with the geometric edge
      itr2 = 0
      for itr2 = 1:mesh.numBC
        if findfirst(mesh.bndry_geo_nums[itr2],geom_face_number) > 0
          break
        end
      end
      start_index = mesh.bndry_offsets[itr2]
      end_index = mesh.bndry_offsets[itr2+1]
      idx_range = start_index:(end_index-1)
      bndry_facenums =sview(mesh.bndryfaces, idx_range) # faces on geometric edge i
      nfaces = length(bndry_facenums)
      for i = 1:nfaces
        bndry_i = bndry_facenums[i]
        # get the local index of the vertices on the boundary face (local face number)
        # vtx_arr = mesh.topo.face_verts[:,bndry_i.face]
        update_coords(mesh, bndry_i.element, real(vertices[itr][:,:,i]))
      end    # End for i = 1:nfaces
    end  # End for itr = 1:length(map.geom_faces)
  end # End

  commit_coords(mesh, sbp, opts)

  return nothing
end
=#

# complexify
import Base.isless

function isless(x::Number, y::Number)

  return isless(real(x), real(y))
end

end # module


text: _is_incremental(::Type{Hypatia.Optimizer}) = false

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: MatthewStuber/EAGO


text:
using ClangCompiler
using ClangCompiler.LLVM

# source file
src = joinpath(@__DIR__, "template.cpp")

# compilation flags
args = get_compiler_args()

# create JIT, generate LLVM IR, and call function
jit = LLJIT(;tm=JITTargetMachine())
irgen = IRGenerator(src, args)
cc = CXCompiler(irgen, jit)
link_process_symbols(cc)
compile(cc)

addr = lookup(jit, "main")
@eval main() = ccall($(pointer(addr)), Cint, ())

main()

# clean up
dispose(cc)


text:
# Compare loading one file in Python vs Julia
using Dates
using InteractiveUtils: versioninfo

if length(ARGS) != 3
    println("Usage: $PROGRAM_FILE <filename> <count> <outputdir>")
	exit()
end

using SASLib, BenchmarkTools
using Humanize: datasize

file = ARGS[1]
cnt = ARGS[2]
dir = ARGS[3]
basename = split(file, "/")[end]
shortname = replace(basename, r"\.sas7bdat" => "")
output = "$dir/py_jl_$(shortname)_$(cnt).md"

prt(msg...) = println(now(), " ", msg...)

# run python part
# result is minimum time in seconds
prt("Running python test for file ", file, " ", cnt, " times")
pyver_cmd = `python -V`
pyres_cmd = `python perf_test1.py $file $cnt`
pyver = readlines(open(pyver_cmd))[1]
pyres = readlines(open(pyres_cmd))[1]  # first line of result is Minimum
py = parse(Float64, match(r"[0-9]*\.[0-9]*", pyres).match)

# read metadata
prt("Reading metadata of data file")
meta = metadata(file)
nrows = meta.nrows
ncols = meta.ncols
nnumcols = count(x -> x == Float64, [ty for (name, ty) in meta.columnsinfo])
nstrcols = ncols - nnumcols

# run julia part
prt("Running julia test")
jlb1 = @benchmark readsas(file, verbose_level=0) samples=parse(Int, cnt)
jl1 = jlb1.times[1] / 1e9   # convert nanoseconds to seconds

# run julia part using regular string array
if nstrcols > 0
    prt("Running julia test (regular string array)")
    jlb2 = @benchmark readsas(file, string_array_fn=Dict(:_all_=>REGULAR_STR_ARRAY),
        verbose_level=0) samples=parse(Int, cnt)
    jl2 = jlb2.times[1] / 1e9   # convert nanoseconds to seconds
    jl = min(jl1, jl2)    # pick faster run
    jltitle = "Julia (ObjectPool)"
else
    jl = jl1
    jltitle = "Julia"
end

# analysis
direction = jl < py ? "faster" : "slower"
ratio = round(direction == "faster" ? py/jl : jl/py, digits = 1)

io = nothing
try
    io = open(output, "w")
    println(io, "# Julia/Python Performance Test Result")
    println(io)
    println(io, "## Summary")
    println(io)
    println(io, "Julia is ~$(ratio)x $direction than Python/Pandas")
    println(io)
    println(io, "## Test File")
    println(io)
    println(io, "Iterations: $cnt")
    println(io)
    println(io, "Filename|Size|Rows|Columns|Numeric Columns|String Columns")
    println(io, "--------|----|----|-------|---------------|--------------")
    println(io, "$basename|$(datasize(lstat(file).size))|$nrows|$ncols|$nnumcols|$nstrcols")
    println(io, "")
    println(io, "## Python")
    println(io, "```")
    println(io, "\$ $(join(pyver_cmd.exec, " "))")
    println(io, pyver)
    println(io, "\$ $(join(pyres_cmd.exec, " "))")
    println(io, pyres)
    println(io, "```")
    println(io)
    println(io, "## $jltitle")
    println(io, "```")
    versioninfo(io)
    println(io)
    show(io, MIME"text/plain"(), jlb1)
    println(io, "\n```")
    if nstrcols > 0
        println(io)
        println(io, "## Julia (Regular String Array)")
        println(io, "```")
        versioninfo(io)
        println(io)
        show(io, MIME"text/plain"(), jlb2)
        println(io, "\n```")
    end
catch err
     println(err)
finally
    io !== nothing && try close(io) catch e println(e) end
end
prt("Written file $output")


text:
"""
A generic interface for two-players, zero-sum games.

Stochastic games and intermediate rewards are supported. By convention,
rewards are expressed from the point of view of the player called _white_.
"""
module GameInterface

export AbstractGame

import ..Util

#####
##### Types
#####

"""
    AbstractGame

Abstract base type for a game environment.

# Constructors

Any subtype `Game` must implement the following constructors:

    Game()

Return an initialized game environment. Note that this constructor does not
have to be deterministic.

    Game(state)

Return a fresh game environment starting at a given state.
"""
abstract type AbstractGame end

function Base.copy(env::Game) where {Game <: AbstractGame}
  return Game(current_state(env))
end

"""
    two_players(::Type{<:AbstractGame}) :: Bool

Return whether or not a game is a two-players game.
"""
function two_players end

"""
    State(Game::Type{<:AbstractGame})

Return the state type corresponding to `Game`.

State objects must be persistent or appear as such as they are stored into
the MCTS tree without copying. They also have to be comparable and hashable.
"""
function State end

"""
    Action(Game::Type{<:AbstractGame})

Return the action type corresponding to `Game`.
"""
function Action end

"""
    actions(::Type{<:AbstractGame})

Return the vector of all game actions.
"""
function actions end

#####
##### Game functions
#####

"""
    white_playing(::Type{<:AbstractGame}, state) :: Bool
    white_playing(env::AbstractGame)
      = white_playing(typeof(env), current_state(env))

Return `true` if white is to play and `false` otherwise. For a one-player
game, it must always return `true`.
"""
function white_playing end

function white_playing(env::AbstractGame)
  return white_playing(typeof(env), current_state(env))
end

"""
    white_reward(env::AbstractGame)

Return the intermediate reward obtained by the white player after the last
transition step. The result is undetermined when called at an initial state.
"""
function white_reward end

"""
    current_state(env::AbstractGame)

Return the game state (which is persistent).
"""
function current_state end

"""
    game_terminated(::AbstractGame)

Return a boolean indicating whether or not the game is in a terminal state.
"""
function game_terminated end

"""
    actions_mask(env::AbstractGame)

Return a boolean mask indicating what actions are available from `env`.

The following identities must hold:

  - `game_terminated(env) || any(actions_mask(env))`
  - `length(actions_mask(env)) == length(actions(typeof(env)))`
"""
function actions_mask end

"""
    play!(env::AbstractGame, action)

Update the game environment by making the current player perform `action`.
Note that this function does not have to be deterministic.
"""
function play! end

"""
    heuristic_value(env::AbstractGame)

Return a heuristic estimate of the state value for the current player.

The given state must be nonfinal and returned values must belong to the
``(-∞, ∞)`` interval.

This function is not needed by AlphaZero but it is useful for building
baselines such as minmax players.
"""
function heuristic_value end

#####
##### Symmetries
#####

"""
    symmetries(::Type{G}, state) where {G <: AbstractGame}

Return the vector of all pairs `(s, σ)` where:
  - `s` is the image of `state` by a nonidentical symmetry
  - `σ` is the associated actions permutation, as an integer vector of
     size `num_actions(Game)`.

A default implementation is provided that returns an empty vector.

# Example

In the game of tic-tac-toe, there are eight symmetries that can be
obtained by composing reflexions and rotations of the board (including the
identity symmetry).
"""
function symmetries(::Type{G}, state) where {G <: AbstractGame}
  return Tuple{State(G), Vector{Int}}[]
end

function test_symmetry(Game, state, (symstate, aperm))
  syms = symmetries
  mask = actions_mask(Game(state))
  symmask = actions_mask(Game(symstate))
  v = falses(length(symmask))
  v[mask] .= true
  v = v[aperm]
  return all(v[symmask]) && !any(v[.~symmask])
end

#####
##### Machine learning interface
#####

"""
    vectorize_state(::Type{<:AbstractGame}, state) :: Array{Float32}

Return a vectorized representation of a state.
"""
function vectorize_state end

#####
##### Interface for interactive exploratory tools
#####

"""
    action_string(::Type{<:AbstractGame}, action) :: String

Return a human-readable string representing the provided action.
"""
function action_string end

"""
    parse_action(::Type{<:AbstractGame}, str::String)

Return the action described by string `str` or `nothing`
if `str` does not denote a valid action.
"""
function parse_action end

"""
    read_state(::Type{G}) where G <: AbstractGame :: Union{State(G), Nothing}

Read a state from the standard input.
Return the corresponding state or `nothing` in case of an invalid input.
"""
function read_state end

"""
    render(env::AbstractGame)

Print the game state on the standard output.
"""
function render end

#####
##### Derived functions
#####

"""
    num_actions(::Type{G})

Return the total number of actions associated with a game.
"""
num_actions(::Type{G}) where G = length(actions(G))

"""
    available_actions(env::AbstractGame)

Return the vector of all available actions.
"""
function available_actions(env::AbstractGame)
  Game = typeof(env)
  mask = actions_mask(env)
  return actions(Game)[mask]
end

"""
    state_dim(::Type{G})

Return a tuple that indicates the shape of a vectorized state representation.
"""
state_dim(::Type{G}) where G = size(vectorize_state(G, current_state(G())))

"""
    apply_random_symmetry(::AbstractGame)

Return a fresh new state that is the image of the given state by a random
symmetry (see [`symmetries`](@ref)).
"""
function apply_random_symmetry(env::Game) where {Game <: AbstractGame}
  symstate, _ = rand(symmetries(Game, current_state(env)))
  return Game(symstate)
end

function state_memsize(::Type{G}) where G
  return Base.summarysize(current_state(G()))
end

end

"""
    GameType(T)

Return the `Game` type associated with an object
(such as a network, a player, an MCTS environment...)
"""
function GameType end


text:
using StaticArrays
using RobotDynamics: KnotPoint, state,control
n,m = 10,5
x = @SVector rand(n)
u = @SVector rand(m)
t = 0.0
dt = 0.1
z = KnotPoint(x,u,t,dt)
RobotDynamics.setstate!(z, 2x)
RobotDynamics.setcontrol!(z, 2u)
@test RD.state(z) ≈ 2x
@test RD.control(z) ≈ 2u
RD.setdata!(z, 3*[x;u])
@test state(z) ≈ 3x
@test control(z) ≈ 3u
@test z isa KnotPoint

# Terminal time step
z = KnotPoint(x,u,t,dt)
@test RD.is_terminal(z) == false
@test RD.getdata(z) ≈ [x;u]
z.dt = 0.0
@test RD.is_terminal(z) == true

# StaticKnotPoint
t = 1.3
z = RD.StaticKnotPoint(x, u, t, dt)
@test_throws ErrorException RobotDynamics.setstate!(z, x)
@test_throws ErrorException RobotDynamics.setcontrol!(z, u)
@test_throws ErrorException RobotDynamics.setdata!(z, [x;u])

z0 = KnotPoint(2x,2u,t,dt)
z = RD.StaticKnotPoint(z0, [x;u])
@test state(z) ≈ x
@test control(z) ≈ u
@test z.dt ≈ dt
@test z.t ≈ 1.3
z = RD.StaticKnotPoint(z0)
@test state(z) ≈ 2x
@test control(z) ≈ 2u


text:
export showinfo

function showinfo(io::IO, a, comment::AbstractString = ""; showheader = true)
    v(a::Number) = [a]
    v(a) = vec(a)
    s(a::AbstractString) = length(a)
    s(a) = size(a)
    med(a::AbstractString) = @p map a uint8 | median | round Int _ | char
    med(a) = median(a)
    if isa(a, Union{Number, AbstractArray, AbstractString})
        showheader && print(io,  length(comment)==0 ? "----  " : comment*"  --  ")
        println(io, "    type: $(typeof(a))   size: $(s(a))")
        try
            if !isa(a, Array) || eltype(a)<:Number
                println(io, "    min:  $(minimum(a))   max: $(maximum(a))\n    mean: $(mean(a))   median: $(med(v(a)))")
            end
        catch
        end
    else
        println(io, "type: $(typeof(a)) ", length(comment)==0 ? "----  " : comment*"  --  ")
    end
    a
end
showinfo(a...; kargs...) = showinfo(stdout, a...; kargs...)


text:
### Day 18
## Like a Rogue
## Author: Thanasis Georgiou
workspace()

# If the neighbors of a tile match any of these then it's trapped
trapScenarios = [
  ['^', '^', '.'],
  ['.', '^', '^'],
  ['^', '.', '.'],
  ['.', '.', '^'],
]

# Check if the neighbors of a tile match any of the trap scenarios
function isTrap(neighbors::Array{Char})::Bool
  for trapScenario in trapScenarios
    if neighbors == trapScenario
      return true
    end
  end
  return false
end

# Puzzle input
input = "......^.^^.....^^^^^^^^^...^.^..^^.^^^..^.^..^.^^^.^^^^..^^.^.^.....^^^^^..^..^^^..^^.^.^..^^..^^^.."
rowLength = length(input)
#rowCount = 40
rowCount = 400000

# Floors plan
floorPlan = Array(Array{Char}, rowCount)
floorPlan[1] = collect(input)
row = 2

# Keep count of how many safe tiles we have found
safeTiles = count(t -> t == '.', floorPlan[1])

# Calculate floor plan
while row <= rowCount
  # Add new row to plan
  floorPlan[row] = Array(Char, rowLength)

  # Calculate first tile
  neighbors = Array(Char, 3)
  neighbors[1] = '.'
  neighbors[2:3] = floorPlan[row - 1][1:2]
  floorPlan[row][1] = isTrap(neighbors)? '^' : '.'
  safeTiles += floorPlan[row][1] == '.'? 1 : 0

  # Calculate last tile
  neighbors[3] = '.'
  neighbors[1:2] = floorPlan[row - 1][(end - 1):end]
  floorPlan[row][end] = isTrap(neighbors)? '^' : '.'
  safeTiles += floorPlan[row][end] == '.'? 1 : 0

  # Calculate the rest
  for i = 2:(rowLength - 1)
    neighbors = floorPlan[row - 1][(i - 1):(i + 1)]
    floorPlan[row][i] = isTrap(neighbors)? '^' : '.'
    safeTiles += floorPlan[row][i] == '.'? 1 : 0
  end

  # Move to next row
  row += 1
end

println("Amount of safe tiles: $safeTiles")


text:
@testset "push integers" begin
    x = IntRangeSet{Int}()
    @test push!(x, 3) == x
    @test collect(x)  == [3:3]
    @test push!(x, 3) == x
    @test collect(x)  == [3:3]
    @test push!(x, 4) == x
    @test collect(x)  == [3:4]
    @test push!(x, 4) == x
    @test collect(x)  == [3:4]
    @test push!(x, 6) == x
    @test collect(x)  == [3:4, 6:6]
    @test push!(x, 5) == x
    @test collect(x)  == [3:6]
    @test push!(x, 0) == x
    @test collect(x)  == [0:0, 3:6]
    @test push!(x, 1) == x
    @test collect(x)  == [0:1, 3:6]
    @test push!(x, 2) == x
    @test collect(x)  == [0:6]
end

@testset "push ranges" begin
    x = IntRangeSet{Int}()
    @test push!(x, 2:5) == x
    @test collect(x)  == [2:5]
    @test push!(x, 15:15) == x
    @test collect(x)  == [2:5, 15:15]
    @test push!(x, 15:14) == x
    @test collect(x)  == [2:5, 15:15]
    @test push!(x, 15:15) == x
    @test collect(x)  == [2:5, 15:15]
    @test push!(x, 1:3) == x
    @test collect(x)  == [1:5, 15:15]
    @test push!(x, 0:0) == x
    @test collect(x)  == [0:5, 15:15]
    @test push!(x, 5:8) == x
    @test collect(x)  == [0:8, 15:15]
    @test push!(x, 14:16) == x
    @test collect(x)  == [0:8, 14:16]
    @test push!(x, 3:5) == x
    @test collect(x)  == [0:8, 14:16]
    @test push!(x, 0:20) == x
    @test collect(x)  == [0:20]
end

@testset "in operator" begin
    x = IntRangeSet{Int}()
    @test (2 in x) == false
    push!(x, 2:4)
    @test (2 in x) == true
    @test (1 in x) == false
    push!(x, 6)
    @test (6 in x) == true
    @test (5 in x) == false
    push!(x, 4:6)
    @test (5 in x) == true
end

@testset "union" begin
    x = IntRangeSet{Int}()
    y = IntRangeSet{Int}()
    push!(x, 2:4)
    push!(y, 6)
    @test collect(union(x,y)) == [2:4, 6:6]
    push!(y, 5:6)
    union!(x, y)
    @test collect(x) == [2:6]
end

@testset "intersect" begin
    x = IntRangeSet{Int}()
    y = IntRangeSet{Int}()
    push!(x, 2:4)
    push!(x, 8:10)
    push!(y, 3:9)
    @test collect(intersect(x, y)) == [3:4, 8:9]
end

@testset "additional random test" begin
    a = IntSet()
    b = IntRangeSet{Int}()

    for i in 1:50_000
        data = rand(1:1_000_000)
        if rand() < .2 # range
            data = data:data+rand(1:5)^2
            foreach(x->push!(a, x), data)
            push!(b, data)
        else # singel integer
            push!(a, data)
            push!(b, data)
        end
    end

    for i in 1:200_000
        data = rand(1:1_000_025)
        @test (data in a) == (data in b)
    end
end


text:
module IrisFlower

using CSV,
	DataFrames,
	Feather,
	Gadfly,
	IJulia,
	JSON,
	Pkg,
	RDatasets,
	Statistics
iris = dataset("datasets", "iris")

end # module


text:
# This file provides getter-type methods which essentially define the fields
# needed in any concrete subtype of SeisModel.

"""
    depth(m, radius) -> depth

Return the `depth` in km for the model `m` given a `radius` in km.
"""
depth(m::SeisModel, radius) = surface_radius(m) - radius

"""
    hasattenuation(m) -> ::Bool

Return `true` if the model `m` includes attenuation and `false` otherwise.
"""
hasattenuation(m::SeisModel) = m.attenuation

"""
	hasdensity(m) -> ::Bool

Return `true` if density is defined for model `m`, and `false` otherwise.
"""
hasdensity(m::SeisModel) = !isempty(m.density)

"""
    hasreffrequency(m::PREMPolyModel) -> ::Bool

Return `true` if a reference frequency is defined for `PREMPolyModel` `m`,
and `false` otherwise.
"""
hasreffrequency(m::PREMPolyModel) = !isnan(m.fref)

"""
    isanisotropic(m) -> ::Bool

Return `true` if the model `m` is anisotropic and `false` otherwise.
"""
isanisotropic(m::SeisModel) = m.aniso

"""
    radius(m, depth) -> radius

Return the `radius` in km for the model `m` given a `depth` in km.
"""
radius(m::SeisModel, depth) = surface_radius(m) - depth

"""
    reffrequency(m::PREMPolyModel) = f

Return the reference frequency for a `PREMPolyModel` `m` in Hz.
"""
function reffrequency(m::PREMPolyModel)
    hasreffrequency(m) ||
        throw(ArgumentError("model does not have a reference frequency"))
    m.fref
end

"""
    surface_radius(m) -> radius

Return the surface radius in km of the `SeisModel` `m`.
"""
surface_radius(m::SeisModel) = m.a

"""
    voigt_kappa_mu(A, C, L, N, F) -> κ, μ

Return the Voigt average bulk modulus `κ` and shear modulus `μ` given the Love parameters
of a transversely isotropic medium.  If A, C, etc., contain density, then κ and μ will be
true moduli, whilst if the Love parameters are provided as velocities, then these
are density-normalised 'moduli' with dimensions of L²⋅T².
"""
voigt_kappa_mu(A, C, L, N, F) = (4A + C + 4F - 4N)/9, (A + C - 2F + 5N + 6L)/15

"""
    voigt_velocities(vpv, vsv, vph, vsh, η) -> vp, vs

Return the Voigt average `vp` and `vs`, in the same units as the input, given four
velocities and the radial anisotropy parameter `η`.
"""
function voigt_velocities(vpv, vsv, vph, vsh, η)
    A = vph^2
    C = vpv^2
    L = vsv^2
    N = vsh^2
    F = η*(A - 2L)
    κ, μ = voigt_kappa_mu(A, C, L, N, F) # In km²/s²
    vp = sqrt(κ + 4/3*μ)
    vs = sqrt(μ)
    vp, vs
end