text:
V471Tau_desc = """
# V471 Tau SPHERE SDI Data

Website: https://github.com/carlgogo/VIP_extras

Author: Carlos Gomez-Gonzalez

## Keys
* `:cube` pre-processed SDI tensor/cube
* `:pa` parallactic angles
* `:wl` wavelengths
* `:psf` off-axis PSF
"""

V471Tau_datadep = DataDep(
    "V471Tau",
    V471Tau_desc,
    [
        "https://rawcdn.githack.com/carlgogo/VIP_extras/a7ec7ac9aeabcb5df73b9ea57270c6d66f9e3f57/datasets/sphere_v471tau_cube.fits",
        "https://rawcdn.githack.com/carlgogo/VIP_extras/a7ec7ac9aeabcb5df73b9ea57270c6d66f9e3f57/datasets/sphere_v471tau_pa.fits",
        "https://rawcdn.githack.com/carlgogo/VIP_extras/a7ec7ac9aeabcb5df73b9ea57270c6d66f9e3f57/datasets/sphere_v471tau_psf.fits",
        "https://rawcdn.githack.com/carlgogo/VIP_extras/a7ec7ac9aeabcb5df73b9ea57270c6d66f9e3f57/datasets/sphere_v471tau_wl.fits"
    ],
    "e6c93b47c401a3747df646abb8311fe9ed3ae128c5a415c8831d6075cf642eee"
)

struct V471TauDataset <: HCIDataset end

@doc V471Tau_desc
const V471Tau = V471TauDataset()

Base.keys(::V471TauDataset) = (:cube, :pa, :psf, :wl)

function Base.getindex(::V471TauDataset, key::Symbol)
    if key === :cube
        getdata(datadep"V471Tau/sphere_v471tau_cube.fits")
    elseif key === :pa
        getdata(datadep"V471Tau/sphere_v471tau_pa.fits")
    elseif key === :psf
        getdata(datadep"V471Tau/sphere_v471tau_psf.fits")
    elseif key === :wl
        getdata(datadep"V471Tau/sphere_v471tau_wl.fits")
    else
        throw(KeyError(key))
    end
end


text:
# import CompScienceMeshes.Simplex

# circumcenter(splx) = circumcenter(splx, dimtype(splx))

# function circumcenter(simplex, dt)

#     F = faces(simplex)
#     c1 = circumcenter(F[1], dimtype(F[1]))
#     c2 = circumcenter(F[2], dimtype(F[2]))

#     n1 = normal(F[1])
#     n2 = normal(F[2])
#     N = [n1 -n2]

#     T = tangents(simplex)

#     c = -c1 + c2
#     s = (T'*N) \ (T'*c)
#     return c1 + n1 * s[1]
# end

# function circumcenter(simplex, ::Type{Val{1}})
#     return (simplex.vertices[1] + simplex.vertices[2])/2
# end

function circumcenter(simplex)

    V = vertices(simplex)
    G = V'*V
    T = coordtype(simplex)
    L = @SMatrix ones(T, size(G,1), 1)
    Z = @SMatrix zeros(T,1,1)

    A = [2G L; L' Z]

    D = vec(sum(V.^2, dims=1))
    z = @SVector zeros(T,1)

    b = [D; z]
    x = A \ b

    c = V*x[1:end-1]

    return c
end

ext: jl


text:
__precompile__

module SMPL

include("smplfuncs.jl");
export createSMPL;
export smpl_lbs;
export viz_smpl;
# export smpl_lbs2;

end # module


text:
function thresh1D(arr, x)
    i=findlast(y -> y<x,arr)
end


text:
export plot3d, plot3d!

function load_polyhedra_mesh()
return quote

using .Polyhedra: Mesh

end end  # quote / function load_polyhedra_mesh()

function load_makie()
return quote

using .Makie: mesh, mesh!
using .Makie: Automatic

end end  # quote / function load_makie()

# helper function for 3D plotting; converts S to a polytope in H-representation
function _plot3d_helper(S::LazySet, backend)
    @assert dim(S) <= 3 "plot3d can only be used to plot sets of dimension three (or lower); " *
        "but the given set is $(dim(S))-dimensional"

    @assert applicable(constraints_list, S) "plot3d requires that the list of constraints of `S`, " *
        "`constraints_list(S)` is applicable; try overapproximating with an `HPolytope` first"

    @assert isbounded(S) "plot3d requires a bounded set"

    P = HPolytope(constraints_list(S))
    remove_redundant_constraints!(P)
    P_poly = polyhedron(P, backend=backend)
    P_poly_mesh = Mesh(P_poly)
    return P_poly_mesh
end

"""
    plot3d(S::LazySet; backend=default_polyhedra_backend(S),
           alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing,
           interpolate=false, linewidth=1, overdraw=false, shading=true,
           transparency=true, visible=true)

Plot a three-dimensional convex set using `Makie`.

### Input

- `S`            -- convex set
- `backend`      -- (optional, default: `default_polyhedra_backend(S)`) polyhedral
                    computations backend
- `alpha`        -- (optional, default: `1.0`) float in `[0,1]`; the alpha or
                    transparency value
- `color`        -- (optional, default: `:blue`) `Symbol` or `Colorant`; the color
                    of the main plot element (markers, lines, etc.) and it can be
                    a color symbol/string like `:red`
- `colormap`     -- (optional, default: `:viridis`) the color map of the main plot;
                    call `available_gradients()` to see what gradients are available,
                    and it can also be used as `[:red, :black]`
- `colorrange`   -- (optional, default: `nothing`, which falls back to
                    `Makie.Automatic()`) a tuple `(min, max)`
                    where `min` and `max` specify the data range to be used for
                    indexing the colormap
- `interpolate`  -- (optional, default: `false`) a bool for heatmap and images,
                    it toggles color interpolation between nearby pixels
- `linewidth`    -- (optional, default: `1`) a number that specifies the width of
                    the line in `line` and `linesegments` plots
- `overdraw`     -- (optional, default: `false`)
- `shading`      -- (optional, default: `true`) a boolean that specifies if shading
                    should be on or not (for meshes)
- `transparency` -- (optional, default: `true`) if `true`, the set is transparent
                    otherwise it is displayed as a solid object
- `visible`      -- (optional, default: `true`) a bool that toggles visibility
                    of the plot

For a complete list of attributes and usage see
[Makie's documentation](http://makie.juliaplots.org/stable/plot-attributes).

### Notes

This plot recipe works by computing the list of constraints of `S` and converting
to a polytope in H-representation. Then, this polytope is transformed with
`Polyhedra.Mesh` and it is plotted using the `mesh` function.

If the function `constraints_list` is not applicable to your set `S`, try
overapproximation first; e.g. via

```julia
julia> Sapprox = overapproximate(S, SphericalDirections(10))

julia> plot3d(Sapprox)
```
The number `10` above corresponds to the number of directions considered; for
better resolution use higher values (but it will take longer).

For efficiency consider using the `CDDLib` backend, as in

```julia
julia> using CDDLib

julia> plot3d(Sapprox, backend=CDDLib.Library())
```

### Examples

The functionality requires *both* `Polyhedra` and `Makie`; so after
loading `LazySets`, do `using Makie, Polyhedra` (or `using Polyhedra, Makie`, the
order doesn't matter).

```julia
julia> using LazySets, Makie, Polyhedra

julia> plot3d(10. * rand(Hyperrectangle, dim=3))

julia> plot3d!(10. * rand(Hyperrectangle, dim=3), color=:red)
```
"""
function plot3d(S::LazySet; backend=default_polyhedra_backend(S),
                alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing, interpolate=false,
                linewidth=1, overdraw=false, shading=true, transparency=true, visible=true)
    require(:Makie; fun_name="plot3d")
    require(:Polyhedra; fun_name="plot3d")

    if colorrange == nothing
        colorrange = Automatic()
    end
    P_poly_mesh = _plot3d_helper(S, backend)
    return mesh(P_poly_mesh, alpha=alpha, color=color, colormap=colormap, colorrange=colorrange,
                interpolate=interpolate, linewidth=linewidth, transparency=transparency, visible=visible)
end

"""
    plot3d!(S::LazySet; backend=default_polyhedra_backend(S),
            alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing, interpolate=false,
            linewidth=1, overdraw=false, shading=true, transparency=true, visible=true)

Plot a three-dimensional convex set using Makie.

### Input

See `plot3d` for the description of the inputs. For a complete list of
attributes and usage see [Makie's
documentation](http://makie.juliaplots.org/stable/plot-attributes).

### Notes

See the documentation of `plot3d` for examples.
"""
function plot3d!(S::LazySet; backend=default_polyhedra_backend(S),
                alpha=1.0, color=:blue, colormap=:viridis, colorrange=nothing, interpolate=false,
                linewidth=1, overdraw=false, shading=true, transparency=true, visible=true)
    require(:Makie; fun_name="plot3d!")
    require(:Polyhedra; fun_name="plot3d!")

    if colorrange == nothing
        colorrange = Automatic()
    end
    P_poly_mesh = _plot3d_helper(S, backend)
    return mesh!(P_poly_mesh, alpha=alpha, color=color, colormap=colormap, colorrange=colorrange,
                 interpolate=interpolate, linewidth=linewidth, transparency=transparency, visible=visible)
end


text:
using PyPlot
PyPlot.ioff()
push!(LOAD_PATH, "/home/berceanu/Development/topo-photon/anc/scripts")
import BP
using PyCall
@pyimport matplotlib.gridspec as gspec

# system parameters
const N = 45
const q = 11
const κ = 0.02
const γ = 0.001
const ν = collect(-3.45:0.001:-0.47)
prm = (1/q,γ,κ);
const r = 11 # points in MBZ
# Nk should be an odd multiple of q
const Nk = r*q # zero-padded system size = no of k points
const l = div(Nk-1,2)
const x = -l:l
const δk = 2π/Nk # resolution in mom space
# for fft, full BZ
const k = x * δk
randpmp(s::Int) = BP.randpmp(N; seed=s)
sprans = BP.Spectrum(ν,randpmp(1234), :symmetric, prm...)

# default exact spectrum, first 100 eigenvalues
exdef = BP.ExactStates(100, :symmetric, N, 1/q, κ)
# selected states for plotting
βs = [9, 20, 30,
      38, 59, 99]
ηs = βs + 1
sω0s = [exdef.νs[state]::Float64 for state in ηs] # 6 frequencies

# matplotlib parameters
matplotlib["rcParams"][:update](Dict("axes.labelsize" => 22,
                                     "axes.titlesize" => 20,
                                     "font.size" => 18,
                                     "legend.fontsize" => 14,
                                     "axes.linewidth" => 1.5,
                                     "font.family" => "serif",
                                     "font.serif" => "Computer Modern Roman",
                                     "xtick.labelsize" => 20,
                                     "xtick.major.size" => 5.5,
                                     "xtick.major.width" => 1.5,
                                     "ytick.labelsize" => 20,
                                     "ytick.major.size" => 5.5,
                                     "ytick.major.width" => 1.5,
                                     "text.usetex" => true,
                                     "figure.autolayout" => true))

fig, axes = plt[:subplots](2,3, figsize=(10,7.3))
for i = 1:3 #loop over columns
    # top row
    ax = axes[1,i]
    img = ax[:imshow](abs2(BP.myfft2(BP.getstate(sprans, sω0s[i]),
       k,k)), origin="upper", ColorMap("gist_heat_r"),
       interpolation="none",
       extent=[-π, π, -π, π])
    ax[:set_xticklabels]([])
    ax[:set_xticks]([-π,0,π])
    ax[:set_yticks]([-π,0,π])
    if i == 1 #leftmost panel
        ax[:set_ylabel](L"$k_y$")
        ax[:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    else
        ax[:set_yticklabels]([])
    end
    #bottom row
    ax = axes[2,i]
    ax[:imshow](abs2(BP.myfft2(BP.getstate(sprans, sω0s[i+3]), k,k)),
       origin="upper", ColorMap("gist_heat_r"), interpolation="none",
       extent=[-π, π, -π, π])
    ax[:set_xlabel](L"$k_x$")
    ax[:set_xticks]([-π,0,π])
    ax[:set_yticks]([-π,0,π])
    ax[:set_xticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    if i == 1 # leftmost panel
        ax[:set_ylabel](L"$k_y$")
        ax[:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    else
        ax[:set_yticklabels]([])
    end
end
fig[:savefig]("../../figures/sym_ring.pdf", transparent=true, pad_inches=0.0, bbox_inches="tight")
plt[:close](fig)

# moving the trap
β = 4 # selected state
ω0 = exdef.νs[β + 1]
# trap positions
pos = [(0.0, 0.0) (2.0, 0.0) (5.5, 0.0) (11.0, 0.0);
       (0.0, 0.0) (2.0, 0.0) (5.5, 0.0) (11.0, 0.0);
       (0.0, 2.0) (0.0, 5.5) (0.0, 11.0) (11.0, 11.0)]
bz = Array(Float64, (3,4,Nk,Nk))
for col = 1:4, row = 1:3
    if row == 1 # landau gauge
        spran = BP.Spectrum([ω0], randpmp(1234), :landau, prm...,
           pos[row, col][1], pos[row, col][2])
    else # symmetric gauge
        spran = BP.Spectrum([ω0], randpmp(1234), :symmetric, prm...,
           pos[row, col][1], pos[row, col][2])
    end
    state = BP.getstate(spran, ω0)
    statefft = BP.myfft2(state, k,k)
    bz[row, col, :, :] = abs2(statefft)
end

fig = plt[:figure]()
gs = gspec.GridSpec(3, 4)
gs[:update](top=0.99, bottom=0.08, left=0.07, right=0.99, hspace=0.05)
axes = Array(PyObject, (3,4))
for col = 0:3, row = 0:2
    axes[row + 1, col + 1] = plt[:subplot](get(gs, (row, col)))
end
for col = 1:4, row = 1:3
    axes[row, col][:imshow](squeeze(bz[row, col, :, :], (1,2)),
                            origin="upper", ColorMap("gist_heat_r"),
                            interpolation="none", extent=[-π, π, -π, π])
    axes[row, col][:set_xticks]([-π,0,π])
    axes[row, col][:set_yticks]([-π,0,π])
end
# disable tick labels for bulk panels
for col = 2:4, row = 1:2
    axes[row, col][:set_yticklabels]([])
    axes[row, col][:set_xticklabels]([])
end
# left margin
for row = 1:2
    axes[row, 1][:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    axes[row, 1][:set_xticklabels]([])
    axes[row, 1][:set_ylabel](L"$k_y$", labelpad=-9)
end
# bottom margin
for col = 2:4
    axes[3, col][:set_xticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
    axes[3, col][:set_yticklabels]([])
    axes[3, col][:set_xlabel](L"$k_x$", labelpad=-3)
end
# bottom left corner
axes[3, 1][:set_xticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
axes[3, 1][:set_yticklabels]([L"$-\pi$",L"$0$",L"$\pi$"])
axes[3, 1][:set_xlabel](L"$k_x$", labelpad=-3)
axes[3, 1][:set_ylabel](L"$k_y$", labelpad=-9)
fig[:savefig]("../../figures/fringe_trap.pdf", transparent=true, pad_inches=0.0, bbox_inches="tight")
plt[:close](fig)

max_stars_repo_name: berceanu/topo-photon

max_issues_repo_name: berceanu/topo-photon

max_forks_repo_name: berceanu/topo-photon


text:
constraints = [TruncateQuantiles(0.1, 0.9)
    for i in 1:length(example_uvals)] #i % 2 == 0 ? TruncateStd(0.3) :

# Draw one realisation
@test resample(example_uvals) isa Vector{<:Real}
@test resample(example_uvals, constraints[1]) isa Vector{<:Real}
@test resample(example_uvals, constraints) isa Vector{<:Real}

# Draw multiple realisations
@test resample(example_uvals, 5) isa Vector{Vector{T}} where T<:Real
@test resample(example_uvals, constraints, 5) isa Vector{Vector{T}} where T<:Real
@test resample(example_uvals, constraints[1], 5) isa Vector{Vector{T}} where T<:Real


text:
using OrdinaryDiffEq, DiffEqDevTools, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

prob = prob_ode_linear
sol =solve(prob,Rosenbrock32())
dt₀ = sol.t[2]

prob = prob_ode_2Dlinear
sol =solve(prob,ExplicitRK(tableau=constructBogakiShampine3()))
dt₀ = sol.t[2]

@test  1e-7 < dt₀ < .1
@test_throws ErrorException sol = solve(prob,Euler())
#dt₀ = sol.t[2]

sol3 =solve(prob,ExplicitRK(tableau=constructDormandPrince8_64bit()))
dt₀ = sol3.t[2]

@test 1e-7 < dt₀ < .3

T = Float32
u0 = T.([1.0;0.0;0.0])

tspan = T.((0,100))
prob = remake(prob, u0=u0, tspan=tspan)
@test_nowarn solve(prob, Euler(); dt=T(0.0001))

tspan = T.((2000,2100))
prob = remake(prob, tspan=tspan)
@test_throws ArgumentError solve(prob, Euler(); dt=T(0.0001)) # Loops forever

function rober(du,u,p,t)
  y₁,y₂,y₃ = u
  k₁,k₂,k₃ = p
  du[1] = -k₁*y₁+k₃*y₂*y₃
  du[2] =  k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
  du[3] =  k₂*y₂^2
  nothing
end
u0 = Float32[1.0,0.0,0.0]
tspan = (0f0, 1f5)
params = (4f-2,3f7,1f4)
prob = ODEProblem(rober,u0,tspan,params)
sol = solve(prob, Rosenbrock23())


text:
# Configuration for nearly neutral simuation
const nn_simtype = 1    # nearly neutral infinite alleles
const type_str = "N3"
const mu_list_flag = true
@everywhere const N_list     = [25,50,100,200,400,800,1600,3200,6400]     # popsize
const N_mu_list= [0.5,1.0,2.0,4.0]   # Mutation rate as a multiple of 1.0/N
const mu_list= [0.025,0.01,0.005,0.0025]   # Mutation rate
const L = 1           # number of loci
const ngens = 10000      # Generations after burn-in
const fix_minimum = 0.45  # Fraction of population needed to "fix"-
const burn_in= 4.0    # generations of burn_in as a multiple of 1/mu
dfe = dfe_neutral
dfe_str = "dfe_neutral"

const use_poplist=false   # if true, save a list of populations as a check of correctness

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
function detect_anagrams(subject, candidates)
    base_letters = sort ∘ collect ∘ lowercase
    is_anagram(word) = lowercase(word) != lowercase(subject) && base_letters(word) == base_letters(subject)
    filter(is_anagram, candidates)
end

max_stars_repo_name: koladilip/julia-exercism

max_issues_repo_name: koladilip/julia-exercism

max_forks_repo_name: koladilip/julia-exercism


text:
"""
    replace_code(s)

"""
function replace_code(s::String)

    lines = split(s,"\r\n")

    incodeblock = false

    out = String[]
    for line in lines
        if strip(line) == "{{{"
            incodeblock=true
            push!(out,"```bash")
        elseif strip(line) == "}}}"
            incodeblock=false
            push!(out,"```")
        elseif incodeblock
            push!(out,line)
        elseif !incodeblock

           line = replace_headers(line)

           # Single line code blocks
           line = replace(line, r"\{\{\{(.*?)\}\}\}" => s"`\1`")

           push!(out,line)
        end
    end

    return join(out,"\r\n")
end


text:
include("../shared_source/boilerplate.jl")

import ClimateMachine.BalanceLaws:
    vars_state,
    init_state_prognostic!,
    init_state_auxiliary!,
    compute_gradient_argument!,
    compute_gradient_flux!,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_conditions,
    boundary_state!
import ClimateMachine.DGMethods: DGModel
import ClimateMachine.NumericalFluxes: numerical_flux_first_order!
import ClimateMachine.Orientations: vertical_unit_vector

"""
    ThreeDimensionalCompressibleNavierStokesEquations <: BalanceLaw
A `BalanceLaw` for shallow water modeling.
write out the equations here
# Usage
    ThreeDimensionalCompressibleNavierStokesEquations()
"""
abstract type AbstractFluid3D <: AbstractFluid end
struct Fluid3D <: AbstractFluid3D end

struct ThreeDimensionalCompressibleNavierStokesEquations{
    I,
    D,
    O,
    A,
    T,
    C,
    F,
    BC,
    FT,
} <: AbstractFluid3D
    initial_value_problem::I
    domain::D
    orientation::O
    advection::A
    turbulence::T
    coriolis::C
    forcing::F
    boundary_conditions::BC
    cₛ::FT
    ρₒ::FT
    function ThreeDimensionalCompressibleNavierStokesEquations{FT}(
        initial_value_problem::I,
        domain::D,
        orientation::O,
        advection::A,
        turbulence::T,
        coriolis::C,
        forcing::F,
        boundary_conditions::BC;
        cₛ = FT(sqrt(10)),  # m/s
        ρₒ = FT(1),  #kg/m³
    ) where {FT <: AbstractFloat, I, D, O, A, T, C, F, BC}
        return new{I, D, O, A, T, C, F, BC, FT}(
            initial_value_problem,
            domain,
            orientation,
            advection,
            turbulence,
            coriolis,
            forcing,
            boundary_conditions,
            cₛ,
            ρₒ,
        )
    end
end

CNSE3D = ThreeDimensionalCompressibleNavierStokesEquations

function vars_state(m::CNSE3D, ::Prognostic, T)
    @vars begin
        ρ::T
        ρu::SVector{3, T}
        ρθ::T
    end
end

function init_state_prognostic!(m::CNSE3D, state::Vars, aux::Vars, localgeo, t)
    cnse_init_state!(m, state, aux, localgeo, t)
end

# default initial state if IVP == nothing
function cnse_init_state!(model::CNSE3D, state, aux, localgeo, t)

    x = aux.x
    y = aux.y
    z = aux.z

    ρ = model.ρₒ
    state.ρ = ρ
    state.ρu = ρ * @SVector [-0, -0, -0]
    state.ρθ = ρ

    return nothing
end

# user defined initial state
function cnse_init_state!(
    model::CNSE3D{<:InitialValueProblem},
    state,
    aux,
    localgeo,
    t,
)
    x = aux.x
    y = aux.y
    z = aux.z

    params = model.initial_value_problem.params
    ic = model.initial_value_problem.initial_conditions

    state.ρ = ic.ρ(params, x, y, z)
    state.ρu = ic.ρu(params, x, y, z)
    state.ρθ = ic.ρθ(params, x, y, z)

    return nothing
end

function vars_state(m::CNSE3D, st::Auxiliary, T)
    @vars begin
        x::T
        y::T
        z::T
        orientation::vars_state(m.orientation, st, T)
    end
end

function init_state_auxiliary!(
    m::CNSE3D,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    # update the geopotential Φ in state_auxiliary.orientation.Φ
    init_state_auxiliary!(
        m,
        (m, aux, tmp, geom) ->
            orientation_nodal_init_aux!(m.orientation, m.domain, aux, geom),
        state_auxiliary,
        grid,
        direction,
    )

    # update ∇Φ in state_auxiliary.orientation.∇Φ
    orientation_gradient(m, m.orientation, state_auxiliary, grid, direction)

    # store coordinates and potentially other stuff
    init_state_auxiliary!(
        m,
        (m, aux, tmp, geom) -> cnse_init_aux!(m, aux, geom),
        state_auxiliary,
        grid,
        direction,
    )

    return nothing
end

function orientation_gradient(
    model::CNSE3D,
    ::Orientation,
    state_auxiliary,
    grid,
    direction,
)
    auxiliary_field_gradient!(
        model,
        state_auxiliary,
        ("orientation.∇Φ",),
        state_auxiliary,
        ("orientation.Φ",),
        grid,
        direction,
    )

    return nothing
end

function orientation_gradient(::CNSE3D, ::NoOrientation, _...)
    return nothing
end

function orientation_nodal_init_aux!(
    ::SphericalOrientation,
    domain::Tuple,
    aux::Vars,
    geom::LocalGeometry,
)
    norm_R = norm(geom.coord)
    @inbounds aux.orientation.Φ = norm_R - domain[1]

    return nothing
end

"""
function orientation_nodal_init_aux!(
    ::SuperSphericalOrientation,
    domain::Tuple,
    aux::Vars,
    geom::LocalGeometry,
)
    norm_R = norm(geom.coord)
    @inbounds aux.orientation.Φ = 1 / norm_R^2
end
"""

function orientation_nodal_init_aux!(
    ::FlatOrientation,
    domain::Tuple,
    aux::Vars,
    geom::LocalGeometry,
)
    @inbounds aux.orientation.Φ = geom.coord[3]

    return nothing
end

function orientation_nodal_init_aux!(
    ::NoOrientation,
    domain::Tuple,
    aux::Vars,
    geom::LocalGeometry,
)
    return nothing
end

function cnse_init_aux!(::CNSE3D, aux, geom)
    @inbounds begin
        aux.x = geom.coord[1]
        aux.y = geom.coord[2]
        aux.z = geom.coord[3]
    end

    return nothing
end

function vars_state(m::CNSE3D, ::Gradient, T)
    @vars begin
        ∇ρ::T
        ∇u::SVector{3, T}
        ∇θ::T
        ∇p::T
    end
end

function compute_gradient_argument!(
    model::CNSE3D,
    grad::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    cₛ = model.cₛ
    ρₒ = model.ρₒ

    grad.∇p = (cₛ * ρ)^2 / (2 * ρₒ)

    compute_gradient_argument!(model.turbulence, grad, state, aux, t)
end

@inline function compute_gradient_argument!(
    ::ConstantViscosity,
    grad::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    u = ρu
    θ = ρθ

    grad.∇ρ = ρ
    grad.∇u = u
    grad.∇θ = θ

    return nothing
end

function vars_state(m::CNSE3D, ::GradientFlux, T)
    @vars begin
        μ∇ρ::SVector{3, T}
        ν∇u::SMatrix{3, 3, T, 9}
        κ∇θ::SVector{3, T}
        ∇p::SVector{3, T}
    end
end

function compute_gradient_flux!(
    model::CNSE3D,
    gradflux::Vars,
    grad::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    gradflux.∇p = grad.∇p

    compute_gradient_flux!(
        model,
        model.turbulence,
        gradflux,
        grad,
        state,
        aux,
        t,
    )
end

@inline function compute_gradient_flux!(
    ::CNSE3D,
    turb::ConstantViscosity,
    gradflux::Vars,
    grad::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    μ = turb.μ * I
    ν = turb.ν * I
    κ = turb.κ * I

    gradflux.μ∇ρ = -μ * grad.∇ρ
    gradflux.ν∇u = -ν * grad.∇u
    gradflux.κ∇θ = -κ * grad.∇θ

    return nothing
end

@inline function flux_first_order!(
    model::CNSE3D,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
    direction,
)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    cₛ = model.cₛ
    ρₒ = model.ρₒ

    flux.ρ += ρu
    # flux.ρu += (cₛ * ρ)^2 / (2 * ρₒ) * I

    advective_flux!(model, model.advection, flux, state, aux, t)

    return nothing
end

advective_flux!(::CNSE3D, ::Nothing, _...) = nothing

@inline function advective_flux!(
    ::CNSE3D,
    ::NonLinearAdvectionTerm,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    ρu = state.ρu
    ρθ = state.ρθ

    flux.ρu += ρu ⊗ ρu / ρ
    flux.ρθ += ρu * ρθ / ρ

    return nothing
end

function flux_second_order!(
    model::CNSE3D,
    flux::Grad,
    state::Vars,
    gradflux::Vars,
    ::Vars,
    aux::Vars,
    t::Real,
)
    flux_second_order!(model, model.turbulence, flux, state, gradflux, aux, t)
end

@inline function flux_second_order!(
    ::CNSE3D,
    ::ConstantViscosity,
    flux::Grad,
    state::Vars,
    gradflux::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρ += gradflux.μ∇ρ
    flux.ρu += gradflux.ν∇u
    flux.ρθ += gradflux.κ∇θ

    return nothing
end

@inline function source!(
    model::CNSE3D,
    source::Vars,
    state::Vars,
    gradflux::Vars,
    aux::Vars,
    t::Real,
    direction,
)

    source.ρu -= gradflux.∇p

    coriolis_force!(model, model.coriolis, source, state, aux, t)
    forcing_term!(model, model.forcing, source, state, aux, t)

    return nothing
end

coriolis_force!(::CNSE3D, ::Nothing, _...) = nothing

@inline function coriolis_force!(
    model::CNSE3D,
    coriolis::fPlaneCoriolis,
    source,
    state,
    aux,
    t,
)
    # f × u
    f = [-0, -0, coriolis_parameter(model, coriolis, aux.coords)]
    ρu = state.ρu

    source.ρu -= f × ρu

    return nothing
end

@inline function coriolis_force!(
    model::CNSE3D,
    coriolis::SphereCoriolis,
    source,
    state,
    aux,
    t,
)
    # f × u
    Ω = coriolis.Ω
    f = @SVector [-0, -0, 2Ω]
    ρu = state.ρu
    source.ρu -= f × ρu
    return nothing
end

forcing_term!(::CNSE3D, ::Nothing, _...) = nothing

@inline function forcing_term!(
    model::CNSE3D,
    buoy::Buoyancy,
    source,
    state,
    aux,
    t,
)
    α = buoy.α
    g = buoy.g
    ρθ = state.ρθ

    # as temperature increase, density decreases
    B = -α * g * ρθ
    k̂ = vertical_unit_vector(model.orientation, aux)

    # gravity points downward
    source.ρu -= B * k̂
end

@inline vertical_unit_vector(::Orientation, aux) = aux.orientation.∇Φ
@inline vertical_unit_vector(::NoOrientation, aux) = @SVector [0, 0, 1]

@inline wavespeed(m::CNSE3D, _...) = m.cₛ

roe_average(ρ⁻, ρ⁺, var⁻, var⁺) =
    (sqrt(ρ⁻) * var⁻ + sqrt(ρ⁺) * var⁺) / (sqrt(ρ⁻) + sqrt(ρ⁺))

function numerical_flux_first_order!(
    ::RoeNumericalFlux,
    model::CNSE3D,
    fluxᵀn::Vars{S},
    n⁻::SVector,
    state⁻::Vars{S},
    aux⁻::Vars{A},
    state⁺::Vars{S},
    aux⁺::Vars{A},
    t,
    direction,
) where {S, A}
    numerical_flux_first_order!(
        CentralNumericalFluxFirstOrder(),
        model,
        fluxᵀn,
        n⁻,
        state⁻,
        aux⁻,
        state⁺,
        aux⁺,
        t,
        direction,
    )

    FT = eltype(fluxᵀn)

    # constants and normal vectors
    cₛ = model.cₛ
    ρₒ = model.ρₒ

    # - states
    ρ⁻ = state⁻.ρ
    ρu⁻ = state⁻.ρu
    ρθ⁻ = state⁻.ρθ

    # constructed states
    u⁻ = ρu⁻ / ρ⁻
    θ⁻ = ρθ⁻ / ρ⁻
    uₙ⁻ = u⁻' * n⁻

    # in general thermodynamics
    p⁻ = (cₛ * ρ⁻)^2 / (2 * ρₒ)
    c⁻ = cₛ * sqrt(ρ⁻ / ρₒ)

    # + states
    ρ⁺ = state⁺.ρ
    ρu⁺ = state⁺.ρu
    ρθ⁺ = state⁺.ρθ

    # constructed states
    u⁺ = ρu⁺ / ρ⁺
    θ⁺ = ρθ⁺ / ρ⁺
    uₙ⁺ = u⁺' * n⁻

    # in general thermodynamics
    p⁺ = (cₛ * ρ⁺)^2 / (2 * ρₒ)
    c⁺ = cₛ * sqrt(ρ⁺ / ρₒ)

    # construct roe averges
    ρ = sqrt(ρ⁻ * ρ⁺)
    u = roe_average(ρ⁻, ρ⁺, u⁻, u⁺)
    θ = roe_average(ρ⁻, ρ⁺, θ⁻, θ⁺)
    c = roe_average(ρ⁻, ρ⁺, c⁻, c⁺)

    # construct normal velocity
    uₙ = u' * n⁻

    # differences
    Δρ = ρ⁺ - ρ⁻
    Δp = p⁺ - p⁻
    Δu = u⁺ - u⁻
    Δρθ = ρθ⁺ - ρθ⁻
    Δuₙ = Δu' * n⁻

    # constructed values
    c⁻² = 1 / c^2
    w1 = abs(uₙ - c) * (Δp - ρ * c * Δuₙ) * 0.5 * c⁻²
    w2 = abs(uₙ + c) * (Δp + ρ * c * Δuₙ) * 0.5 * c⁻²
    w3 = abs(uₙ) * (Δρ - Δp * c⁻²)
    w4 = abs(uₙ) * ρ
    w5 = abs(uₙ) * (Δρθ - θ * Δp * c⁻²)

    # fluxes!!!
    fluxᵀn.ρ -= (w1 + w2 + w3) * 0.5
    fluxᵀn.ρu -=
        (
            w1 * (u - c * n⁻) +
            w2 * (u + c * n⁻) +
            w3 * u +
            w4 * (Δu - Δuₙ * n⁻)
        ) * 0.5
    fluxᵀn.ρθ -= ((w1 + w2) * θ + w5) * 0.5

    return nothing
end

boundary_conditions(model::CNSE3D) = model.boundary_conditions

"""
    boundary_state!(nf, ::CNSE3D, args...)
applies boundary conditions for the hyperbolic fluxes
dispatches to a function in CNSEBoundaryConditions
"""
@inline function boundary_state!(nf, bc, model::CNSE3D, args...)
    return _cnse_boundary_state!(nf, bc, model, args...)
end

"""
    cnse_boundary_state!(nf, bc::FluidBC, ::CNSE3D)
splits boundary condition application into velocity
"""
@inline function cnse_boundary_state!(nf, bc::FluidBC, m::CNSE3D, args...)
    cnse_boundary_state!(nf, bc.momentum, m, m.turbulence, args...)
    cnse_boundary_state!(nf, bc.temperature, m, args...)

    return nothing
end

include("bc_momentum.jl")
include("bc_temperature.jl")

"""
STUFF FOR ANDRE'S WRAPPERS
"""

function get_boundary_conditions(
    model::SpatialModel{BL},
) where {BL <: AbstractFluid3D}
    bcs = model.boundary_conditions

    west_east = (check_bc(bcs, :west), check_bc(bcs, :east))
    south_north = (check_bc(bcs, :south), check_bc(bcs, :north))
    bottom_top = (check_bc(bcs, :bottom), check_bc(bcs, :top))

    return (west_east..., south_north..., bottom_top...)
end

function DGModel(
    model::SpatialModel{BL};
    initial_conditions = nothing,
) where {BL <: AbstractFluid3D}
    params = model.parameters
    physics = model.physics

    Lˣ, Lʸ, Lᶻ = length(model.grid.domain)
    bcs = get_boundary_conditions(model)
    FT = eltype(model.grid.numerical.vgeo)

    if !isnothing(initial_conditions)
        initial_conditions = InitialValueProblem(params, initial_conditions)
    end

    balance_law = CNSE3D{FT}(
        initial_conditions,
        (Lˣ, Lʸ, Lᶻ),
        physics.orientation,
        physics.advection,
        physics.dissipation,
        physics.coriolis,
        physics.buoyancy,
        bcs,
        ρₒ = params.ρₒ,
        cₛ = params.cₛ,
    )

    numerical_flux_first_order = model.numerics.flux # should be a function

    rhs = DGModel(
        balance_law,
        model.grid.numerical,
        numerical_flux_first_order,
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    return rhs
end


text:
abstract type ReproducingKernel end
abstract type ReproducingKernel_0 <: ReproducingKernel end
abstract type ReproducingKernel_1 <: ReproducingKernel_0 end
abstract type ReproducingKernel_2 <: ReproducingKernel_1 end

@doc raw"
`struct RK_H0{T} <: ReproducingKernel_0`

Defines a type of reproducing kernel of Bessel Potential space ``H^{n/2 + 1/2}_ε (R^n)`` ('Basic Matérn kernel'):
```math
V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|) \, .
```
# Fields
- `ε::T`: 'scaling parameter' from the Bessel Potential space definition,
           it may be omitted in the struct constructor otherwise it must be greater than zero
"
struct RK_H0{T} <: ReproducingKernel_0
    ε::T
    RK_H0() = new{Float64}(0.0)
    RK_H0(ε) = (@assert ε > 0; new{typeof(float(ε))}(float(ε)))
    RK_H0{T}(ε) where {T} = (@assert ε > 0; new{T}(T(ε)))
end
Base.eltype(::RK_H0{T}) where {T} = T

@doc raw"
`struct RK_H1{T} <: ReproducingKernel_1`

Defines a type of reproducing kernel of Bessel Potential space ``H^{n/2 + 3/2}_ε (R^n)`` ('Linear Matérn kernel'):
```math
V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|)
             (1 + \varepsilon |\xi  - \eta|) \, .
```
# Fields
- `ε::T`: 'scaling parameter' from the Bessel Potential space definition,
           it may be omitted in the struct constructor otherwise it must be greater than zero
"
struct RK_H1{T} <: ReproducingKernel_1
    ε::T
    RK_H1() = new{Float64}(0.0)
    RK_H1(ε) = (@assert ε > 0; new{typeof(float(ε))}(float(ε)))
    RK_H1{T}(ε) where {T} = (@assert ε > 0; new{T}(T(ε)))
end
Base.eltype(::RK_H1{T}) where {T} = T

@doc raw"
`struct RK_H2{T} <: ReproducingKernel_2`

Defines a type of reproducing kernel of Bessel Potential space ``H^{n/2 + 5/2}_ε (R^n)`` ('Quadratic Matérn kernel'):
```math
V(\eta , \xi, \varepsilon) = \exp (-\varepsilon |\xi - \eta|)
             (3 + 3\varepsilon |\xi  - \eta| + \varepsilon ^2 |\xi - \eta| ^2 ) \, .
```
# Fields
- `ε::T`: 'scaling parameter' from the Bessel Potential space definition,
           it may be omitted in the struct constructor otherwise it must be greater than zero
"
struct RK_H2{T} <: ReproducingKernel_2
    ε::T
    RK_H2() = new{Float64}(0.0)
    RK_H2(ε) = (@assert ε > 0; new{typeof(float(ε))}(float(ε)))
    RK_H2{T}(ε) where {T} = (@assert ε > 0; new{T}(T(ε)))
end
Base.eltype(::RK_H2{T}) where {T} = T

@inline function _rk(kernel::RK_H2, η::SVector, ξ::SVector)
    x = kernel.ε * norm(η - ξ)
    return (3 + x * (3 + x)) * exp(-x)
end

@inline function _rk(kernel::RK_H1, η::SVector, ξ::SVector)
    x = kernel.ε * norm(η - ξ)
    return (1 + x) * exp(-x)
end

@inline function _rk(kernel::RK_H0, η::SVector, ξ::SVector)
    x = kernel.ε * norm(η - ξ)
    return exp(-x)
end

@inline function _∂rk_∂e(kernel::RK_H2, η::SVector, ξ::SVector, e::SVector)
    t = η - ξ
    x = kernel.ε * norm(t)
    return kernel.ε^2 * exp(-x) * (1 + x) * (t ⋅ e)
end

@inline function _∂rk_∂e(kernel::RK_H1, η::SVector, ξ::SVector, e::SVector)
    t = η - ξ
    x = kernel.ε * norm(t)
    return kernel.ε^2 * exp(-x) * (t ⋅ e)
end

@inline function _∂rk_∂η(kernel::RK_H2, η::SVector, ξ::SVector)
    t = η - ξ
    x = kernel.ε * norm(t)
    return -kernel.ε^2 * exp(-x) * (1 + x) * t
end

@inline function _∂rk_∂η(kernel::RK_H1, η::SVector, ξ::SVector)
    t = η - ξ
    x = kernel.ε * norm(t)
    return -kernel.ε^2 * exp(-x) * t
end

@inline function _∂rk_∂η(kernel::RK_H0, η::SVector, ξ::SVector)
    # Note: Derivative of spline built with reproducing kernel RK_H0 does not exist at the spline nodes
    t    = η - ξ
    tnrm = norm(t)
    x    = kernel.ε * tnrm
    ∇    = kernel.ε * exp(-x)
    t    = ifelse(x > eps(typeof(x)), t, zeros(t))
    ∇   *= ifelse(x > eps(typeof(x)), inv(tnrm), one(tnrm))
    ∇   *= t
end

@inline function _∂²rk_∂²e(kernel::RK_H2, η::SVector{n}, ξ::SVector{n}, êη::SVector{n}, êξ::SVector{n}) where {n}
    ε     = kernel.ε
    ε²    = ε * ε
    t     = η - ξ
    tnrm  = norm(t)
    x     = ε * tnrm
    ε²e⁻ˣ = ε² * exp(-x)
    ∇²    = ((1 + x) * ε²e⁻ˣ) * (êξ ⋅ êη)
    ∇²   -= (ε² * ε²e⁻ˣ) * (êξ ⋅ t) * (t ⋅ êη)
end

@inline function _∂²rk_∂²e(kernel::RK_H1, η::SVector{n}, ξ::SVector{n}, êη::SVector{n}, êξ::SVector{n}) where {n}
    # Note: Second derivative of spline built with reproducing kernel RK_H1 does not exist at the spline nodes
    ε     = kernel.ε
    ε²    = ε * ε
    t     = η - ξ
    tnrm  = norm(t)
    x     = ε * tnrm
    ε²e⁻ˣ = ε² * exp(-x)
    ∇²    = ε²e⁻ˣ * (êξ ⋅ êη)
    ∇²   -= ifelse(x > eps(typeof(x)), (ε * ε²e⁻ˣ / tnrm) * (êξ ⋅ t) * (t ⋅ êη), zero(∇²))
end

@inline function _∂²rk_∂η∂ξ(kernel::RK_H2, η::SVector{n}, ξ::SVector{n}) where {n}
    ε     = kernel.ε
    ε²    = ε * ε
    t     = η - ξ
    tnrm  = norm(t)
    x     = ε * tnrm
    ε²e⁻ˣ = ε² * exp(-x)
    S     = SMatrix{n,n,typeof(x)}
    ∇²    = S(((1 + x) * ε²e⁻ˣ) * I)
    ∇²   -= ((ε² * ε²e⁻ˣ) * t) * t'
end

@inline function _∂²rk_∂η∂ξ(kernel::RK_H1, η::SVector{n}, ξ::SVector{n}) where {n}
    # Note: Second derivative of spline built with reproducing kernel RK_H1 does not exist at the spline nodes
    ε     = kernel.ε
    ε²    = ε * ε
    t     = η - ξ
    tnrm  = norm(t)
    x     = ε * tnrm
    ε²e⁻ˣ = ε² * exp(-x)
    S     = SMatrix{n,n,typeof(x)}
    ∇²    = S(ε²e⁻ˣ * I)
    ∇²   -= ifelse(x > eps(typeof(x)), ((ε * ε²e⁻ˣ / tnrm) * t) * t', zeros(S))
end


text:
# basics for keyword system

KParam = Dict{Symbol,Any}

# KPAtom
mutable struct KPAtom
 Psym :: Symbol
 args :: KParam
end

# KPLiteral
mutable struct KPLiteral
 sign :: Any
 atom :: KPAtom
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
#: ==========================================
#: the stupid part of the SmartParser
#: ==========================================

#! DO NOT CHANGE THE SEQUENCE OF THE LIST BELOW
#! IT IS VERYY DELLICATTE

global const QE_KEYWORDS = [
    #: note the white space !!!
    r"Dynamical RAM for\s+(wfc|wfc\(w\.\s+buffer\)|str\.\s*fact|local\s+pot|nlocal\s+pot|qrad|rho\,v\,vnew|rhoin|rho\*nmix|G\-vectors|h\,s\,v\(r\/c\)|\<psi\|beta\>|psi|hpsi|spsi|wfcinit\/wfcrot|addusdens|addusforce|addusstress)" => " __QEDynRAMfor__ ",
    r"(kinetic|local|nonloc\.|hartree|exc\-cor|corecor|ewald|hubbard|london|DFT\-D3|XDM|dft\-nl|TS\-vdW)\s+stress" => " __QEstressKW__ ",
    r"(local|non\-local|ionic|(core|SCF)\s+corr(ection)?|Hubbard|) contrib(\.|ution)\s+to forces" => " __QEforceKW__ ",
    r"(PWSCF|electrons|forces|stress|sum\_band(\:bec)*|c\_bands|init\_run|update\_pot|wfcinit(\:atom|\:wfcr)*|potinit|hinit0|v\_of\_rho|v\_h|v\_xc|newd|mix\_rho|init\_us\_2|addusdens|\*egterg|[csgh]\_psi(\:pot|\:calbec)*|cdiaghg(\:chol|\:inve|\:para)*|cegterg(\:over|\:upda|\:last)*|vloc\_psi|add\_vuspsi|PAW(\_pot|\_symme)*)\s*\:" => " __QEelROUTINES__ ",
    r"(PHONON|phq(_setup|_init)?|init(_vloc|_us_1)?|dynmat(\d|_us)|phqscf|dynmatrix|solve_linter|drhodv|d2ionq|phqscf|dvqpsi_us|ortho|incdrhoscf|vpsifft|dv_of_drho|mix_pot|ef_shift|localdos|psymdvscf|dvqpsi_us(_on)?|cgsolve|last|add_vuspsi|incdrhoscf)\s*\:" => " __QEphROUTINES__ ",
    r"(calbec|fft(s|w|\_scatt\_xy|\_scatt\_yz)?|davcio|write\_rec)\s*\:" => " __QEgenROUTINES__ ",
    r"stopping\s*\.+"                                                    => " __QESTOPPING__ ",

    # really stupid but they are NOT chemical formulae
    [Regex("(?<![0-9A-Za-z_\\-\\.\\/])$i(?=[^0-9A-Za-z_\\.\\/\\-\\n])")  =>  " __KW$(replace(i,r"(\\\.|\\\-|\d)"=>""))__ "
            for i ∈ ["PW","PWSCF","PH","PHONON","UPF","SCF","CPU","PBC","INF","FHI98PP", "pw\\.x", "fhi2upf\\.x", "ld1\\.x", "calypso\\.x", "vasp\\.x"] ] ...
]

global const LAMMPS_KEYWORDS = [
    r"\$\{\w+\}" => " __LAMMPSVAR__ ",
]

##

global const __NUM_TO_STERLING__ = [
    r"\d+(\.\d+)\%" => " __PERCENTAGE__ ",

    #: this substitution fucks up a lot of lines
    #: still in development
    r"(?<![\w])([+-]?\d+(\.\d*)?|[+-]?\d*\.\d+)((e|E)[+-]?\d+)?(?=[^\w]|$)"  => "£",

    #: aftermath of the number substitution
    "(£"   =>  "( £",
    "£)"   =>  "£ )",
    "£|"   =>  "£ |",
    "|£"   =>  "| £",
    "£,"   =>  "£ ,",
    ",£"   =>  ", £",

    "£-£"  =>  "£ - £",
    "£+£"  =>  "£ + £",

    "£-"   =>  "£ - ",
    "-£"   =>  " - £",
    "£/"   =>  "£ / ",
    "/£"   =>  " / £",
    "£*"   =>  "£ * ",
    "*£"   =>  " * £",
    "£:"   =>  " : £",
    ":£"   =>  "£ : ",
]

global const __SINGLE_ELEM__str__ = "(A[cglmrstu]|B[aehikr]?|C[adeflmnorsu]?|D[bsy]|E[rsu]|F[elmr]?|G[ade]|H[efgos]?|I[nr]?|Kr?|L[airuv]|M[dgnot]|N[abdeiop]?|Os?|P[abdmortu]?|R[abefghnu]|S[bcegimnr]?|T[abcehilm]|U(u[opst])?|V|W|Xe|Yb?|Z[nr])"

##//global const __CHEMELEM__rstr__ = Regex("(?<![0-9A-Za-z])$(__SINGLE_ELEM__str__)(?=[^0-9A-Za-z])")

global const __CHEMFORMULA__rstr__  = Regex("(?<![0-9A-Za-z_\\.\\-\\/])($(__SINGLE_ELEM__str__)\\d*)+(?=[^0-9A-Za-z_\\.\\-\\/\\n])")

##

global const  MASK_RULES = [

    QE_KEYWORDS...,

    LAMMPS_KEYWORDS...,

    #: note the white space !!!
    r"MD5\s*check sum\s*:\s*[a-f0-9]{32}"                   => " __CHKSUM__ ",
    r"point\s+group\s+[A-Za-z0-9\_]+(\s*\(m-3m\))?"         => " __POINTGROUP__ ",

    r"Ry\/bohr\*\*3"                                        => " __UNITSTRESS__ ",
    r"Ry\/Bohr"                                             => " __UNITFORCEa__ ",
    r"Ry\/au"                                               => " __UNITFORCEb__ ",
    r"\d\s*pi\/alat"                                        => " __UNITTWOPIALAT__ ",
    r"g\/cm\^3"                                             => " __UNITDENSITY__ ",
    r"(?<![A-Za-z_\-])Ry"                                   => " __Ry__ ",
    r"(?<![\w\d\-\/])\(\s*alat\s*\)(?=[^\w\d\/]|$)"         => " __ALAT__ ",
    r"(?<![\w\d\-\/])ev(?=[^\w\d\/]|$)"                     => " __EV__ ",
    r"(?<![\w\d\-\/])\(bohr\)(?=[^\w\d\/]|$)"               => " __BOHR__ ",

    r"http\:\/(\/[^\^\/\s]+)+(\/)?"                         => " __URL__ ",
    r"(?<![0-9A-Za-z\*])((\/\w[^\^\/\s\$]*)+(\/)?|(\w[^\^\/\s\$]*\/)+)(?=[^\.0-9A-Za-z\(\)]|$)" => " __FULLPATH__ ",

    r"\(a\.u\.\)\^3"                                        => " __UNITVOLa__ ",
    r"a\.u\.\^3"                                            => " __UNITVOLb__ ",
    r"a\.u\."                                               => " __au__ ",

    r"v(\.\d+){2,3}"                                        => " __VERSIONa__ ",
    r"((\d+h\s*)?\d+m\s*)?\d+\.\d+s"                        => " __DURATION__ ",

    #r"[^0-9\^\/\s\.\(\)][^\^\/\s\.\(\)]{1,20}\.([^\^\/\s\.\(\)]{1,20}\.)*[A-Za-z0-9]+"  => " __ABSPATH__ ",
    #r"[^\^\/\s\.\(\)]{1,20}\w\.([^\^\/\s\.\(\)]{1,20}\.)*[A-Za-z][A-Za-z0-9]*"  => " __ABSPATH__ ",
    r"(?<![0-9A-Za-z\*\.])[^\^\/\s\.\(\)]{1,20}\w\.([^\^\/\s\.\(\)]{1,20}\.)*[A-Za-z][A-Za-z0-9]*(?=[^\.0-9A-Za-z\/]|$)" => " __RELPATH__ ",

    r"Ang\^3"                                               => " __UNITVOLc__ ",
    r"kbar"                                                 => " __UNITkbar__ ",
    r"\[cm-1\]"                                             => " __UNITCMINV__ ",

    r"([01]?\d|2[0-3]):([0-5 ]?\d):([0-5 ]?\d)"             => " __HHMMSS__ ",
    r"[1-3 ]\d[A-Za-z]{2,9}(19|20)\d\d"                     => " __DATEa__ ",

    r" P\s*\= "                                             => " __PRESSUREEQS__ ",

    r"\[\s*\-?\d+\s*,\s*\-?\d+\s*,\s*\-?\d+\s*\]"           => " __MILLER__ ",
    r"\(\s*\d+\s*,\s*\d+\s*,\s*\d+\s*\)"                    => " __THREETUPLES__ ",

    __CHEMFORMULA__rstr__ => " __CHEM__ ",
    ##//__CHEMELEM__rstr__ => " __CHEMELEM__ ",

    r"(?<![0-9A-Za-z_\-])[ABTE]\'*\_[1-5]?[ug](?=[^0-9A-Za-z_\-]|$)"                    => " __REPSYMBOL__ " ,
    r"(?<![0-9A-Za-z_\-])([23468]?[CS][2346]'?|i|E|[2346]s_[hd])(?=[^0-9A-Za-z_\-]|$)"  => " __GRPSYMBOL__ " ,
    r"(?<![0-9A-Za-z_\-\/])[1-6][SsPpDdFf](?=[^0-9A-Za-z_\-\/]|$)"                      => " __ATOMORBIT__ ",

    __NUM_TO_STERLING__...,

    r"[A-Za-z][^\_\\\/\%\s]*(\_[A-Za-z][^\_\\\/\%\s]*)+"    => " __SYMBOLtypeA__ ",
    r"[A-Za-z]+\d+[A-Za-z]*"                                => " __SYMBOLtypeB__ ",  #: this is rather unfortunate
    r"[A-Za-z]*\d+[A-Za-z]+"                                => " __SYMBOLtypeC__ ",  #: this is rather unfortunate
]

# these are used in parse_file for performance optimization
global const  MASK_RULES_DIC_INV = Dict(v=>k for (k,v) ∈ MASK_RULES if (k isa Regex))
global const  MASK_RULES_DIC_INV1 = Dict(strip(v)=>k for (k,v) ∈ MASK_RULES if (k isa Regex))
global const  MASK_RULES_DIC_INV_KEYS_STRIPPED_NO_£ = Set([strip.(x) for x in keys(MASK_RULES_DIC_INV) if x!="£"])

global const __preproc__ = [
    r"(\d+)q-points"          => s"\1 q-points",
    r"\#(\s*\d+): "           => s"# \1 : ",
    r"(\w)\)\:"               =>s"\1 ) :",
    r" (\d+\s*)\*(\s*\d+) "   => s" \1 * \2 ",
    r"([\<\>\=\+\-\*])(\d+)([\*\+\-])" => s"\1 \2 \3",
    r"(\d)ns(?=[^0-9A-Za-z])" => s"\1 ns ",
    r"N(\s+)xq"               => s"NUM\1 xq",  # QE PH calculation "N       xq(1) ..."
    r"\(\s*ev\s*\)"           => " ( ev ) ",
    "G-vectors)"              => "G-vectors )",
    #r" hinit(\d+) " => s" hinit \1 ",
    #r" dynmat(\d+) " => s" dynmat \1 ",
    #r" d(\d+)ionq " => s" d \1 ionq ",
    #"h,s,v(r/c)" => "h,s,v( r / c )",
    #"wfcinit/wfcrot" => "wfcinit / wfcrot",
    #"atoms/cell" => "atoms / cell",
    #"proc/nbgrp/npool/nimage" => " proc / nbgrp / npool / nimage ",
]

max_stars_repo_name: algorithmx/SmartParser

max_issues_repo_name: algorithmx/SmartParser

max_forks_repo_name: algorithmx/SmartParser


text:
export l2Norm_primaldual

"""
Force the rows of L and R to have norm at most B
"""
function l2Norm_primaldual(x::Tx, weights, params::Dict{String,<:Any}) where {ETx<:Number, Tx<:AbstractVector{ETx}}

p = norm(x[:])

return p

end


text:
CholType{T,S<:AbstractMatrix} = Cholesky{T,S}
chol_lower(a::Matrix) = cholesky(a).L

if HAVE_CHOLMOD
    CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}

    chol_lower(cf::CholTypeSparse) = cf.L
end


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using JLD2
import Base: flush

mutable struct FileInstance <: Instance
    py::Union{Nothing,PyCall.PyObject}
    loaded::Union{Nothing,JuMPInstance}
    filename::AbstractString
    sample::PyCall.PyObject
    build_model::Function
    mode::String

    function FileInstance(
        filename::AbstractString,
        build_model::Function;
        mode::String = "a",
    )::FileInstance
        instance = new(nothing, nothing, filename, nothing, build_model, mode)
        instance.py = PyFileInstance(instance)
        if mode != "r" || isfile("$filename.h5")
            instance.sample = Hdf5Sample("$filename.h5", mode = mode)
        end
        instance.filename = filename
        return instance
    end
end

function _load!(instance::FileInstance)
    if instance.loaded === nothing
        data = load_data(instance.filename)
        instance.loaded = JuMPInstance(instance.build_model(data))
    end
end

function free(instance::FileInstance)
    instance.loaded = nothing
end

function to_model(instance::FileInstance)
    _load!(instance)
    return to_model(instance.loaded)
end

function get_instance_features(instance::FileInstance)
    _load!(instance)
    return get_instance_features(instance.loaded)
end

function get_variable_features(instance::FileInstance, names)
    _load!(instance)
    return get_variable_features(instance.loaded, names)
end

function get_variable_categories(instance::FileInstance, names)
    _load!(instance)
    return get_variable_categories(instance.loaded, names)
end

function get_constraint_features(instance::FileInstance, names)
    _load!(instance)
    return get_constraint_features(instance.loaded, names)
end

function get_constraint_categories(instance::FileInstance, names)
    _load!(instance)
    return get_constraint_categories(instance.loaded, names)
end

function find_violated_lazy_constraints(instance::FileInstance, solver)
    _load!(instance)
    return find_violated_lazy_constraints(instance.loaded, solver)
end

function enforce_lazy_constraint(instance::FileInstance, solver, violation)
    _load!(instance)
    return enforce_lazy_constraint(instance.loaded, solver, violation)
end

function get_samples(instance::FileInstance)
    return [instance.sample]
end

function create_sample!(instance::FileInstance)
    if instance.mode == "r"
        return MemorySample()
    else
        return instance.sample
    end
end

function save_data(filename::AbstractString, data)::Nothing
    jldsave(filename, data = data)
end

function load_data(filename::AbstractString)
    jldopen(filename, "r") do file
        return file["data"]
    end
end

function load(filename::AbstractString, build_model::Function)
    jldopen(filename, "r") do file
        return build_model(file["data"])
    end
end

function save(data::AbstractVector, dirname::String)::Nothing
    mkpath(dirname)
    for (i, d) in enumerate(data)
        filename = joinpath(dirname, @sprintf("%06d.jld2", i))
        jldsave(filename, data = d)
    end
end

function solve!(
    solver::LearningSolver,
    filenames::Vector,
    build_model::Function;
    tee::Bool = false,
)
    for filename in filenames
        solve!(solver, filename, build_model; tee)
    end
end

function fit!(
    solver::LearningSolver,
    filenames::Vector,
    build_model::Function;
    tee::Bool = false,
)
    instances = [FileInstance(f, build_model) for f in filenames]
    fit!(solver, instances)
end

function solve!(
    solver::LearningSolver,
    filename::AbstractString,
    build_model::Function;
    tee::Bool = false,
)
    solve!(solver, FileInstance(filename, build_model); tee)
end

function __init_PyFileInstance__()
    @pydef mutable struct Class <: miplearn.Instance
        function __init__(self, jl)
            self.jl = jl
        end
        to_model(self) = to_model(self.jl)
        get_instance_features(self) = get_instance_features(self.jl)
        get_variable_features(self, names) =
            get_variable_features(self.jl, from_str_array(names))
        get_variable_categories(self, names) =
            to_str_array(get_variable_categories(self.jl, from_str_array(names)))
        get_samples(self) = get_samples(self.jl)
        create_sample(self) = create_sample!(self.jl)
        find_violated_lazy_constraints(self, solver, _) =
            find_violated_lazy_constraints(self.jl, solver)
        enforce_lazy_constraint(self, solver, _, violation) =
            enforce_lazy_constraint(self.jl, solver, violation)
        free(self) = free(self.jl)

        # FIXME: The two functions below are disabled because they break lazy loading
        # of FileInstance.

        # get_constraint_features(self, names) =
        #     get_constraint_features(self.jl, from_str_array(names))
        # get_constraint_categories(self, names) =
        #     to_str_array(get_constraint_categories(self.jl, from_str_array(names)))

    end
    copy!(PyFileInstance, Class)
end

export FileInstance


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using Geodesy

function calculate_distance(source_lat, source_lon, dest_lat, dest_lon)::Float64
    x = LLA(source_lat, source_lon, 0.0)
    y = LLA(dest_lat, dest_lon, 0.0)
    return round(distance(x, y) / 1000.0, digits = 2)
end

function build_graph(instance::Instance)::Graph
    arcs = []
    next_index = 0
    process_nodes = ProcessNode[]
    plant_shipping_nodes = ShippingNode[]
    collection_shipping_nodes = ShippingNode[]

    name_to_process_node_map = Dict{Tuple{AbstractString,AbstractString},ProcessNode}()

    process_nodes_by_input_product =
        Dict(product => ProcessNode[] for product in instance.products)
    shipping_nodes_by_plant = Dict(plant => [] for plant in instance.plants)

    # Build collection center shipping nodes
    for center in instance.collection_centers
        node = ShippingNode(next_index, center, center.product, [], [])
        next_index += 1
        push!(collection_shipping_nodes, node)
    end

    # Build process and shipping nodes for plants
    for plant in instance.plants
        pn = ProcessNode(next_index, plant, [], [])
        next_index += 1
        push!(process_nodes, pn)
        push!(process_nodes_by_input_product[plant.input], pn)

        name_to_process_node_map[(plant.plant_name, plant.location_name)] = pn

        for product in keys(plant.output)
            sn = ShippingNode(next_index, plant, product, [], [])
            next_index += 1
            push!(plant_shipping_nodes, sn)
            push!(shipping_nodes_by_plant[plant], sn)
        end
    end

    # Build arcs from collection centers to plants, and from one plant to another
    for source in [collection_shipping_nodes; plant_shipping_nodes]
        for dest in process_nodes_by_input_product[source.product]
            distance = calculate_distance(
                source.location.latitude,
                source.location.longitude,
                dest.location.latitude,
                dest.location.longitude,
            )
            values = Dict("distance" => distance)
            arc = Arc(source, dest, values)
            push!(source.outgoing_arcs, arc)
            push!(dest.incoming_arcs, arc)
            push!(arcs, arc)
        end
    end

    # Build arcs from process nodes to shipping nodes within a plant
    for source in process_nodes
        plant = source.location
        for dest in shipping_nodes_by_plant[plant]
            weight = plant.output[dest.product]
            values = Dict("weight" => weight)
            arc = Arc(source, dest, values)
            push!(source.outgoing_arcs, arc)
            push!(dest.incoming_arcs, arc)
            push!(arcs, arc)
        end
    end

    return Graph(
        process_nodes,
        plant_shipping_nodes,
        collection_shipping_nodes,
        arcs,
        name_to_process_node_map,
    )
end

function print_graph_stats(instance::Instance, graph::Graph)::Nothing
    @info @sprintf("    %12d time periods", instance.time)
    @info @sprintf("    %12d process nodes", length(graph.process_nodes))
    @info @sprintf("    %12d shipping nodes (plant)", length(graph.plant_shipping_nodes))
    @info @sprintf(
        "    %12d shipping nodes (collection)",
        length(graph.collection_shipping_nodes)
    )
    @info @sprintf("    %12d arcs", length(graph.arcs))
    return
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


text:
using PyCallChainRules.Torch: TorchModuleWrapper, torch, functorch, ispysetup, pyfrom_dlpack

using Test
using Zygote
using Flux
using ChainRulesCore: NoTangent, AbstractZero
import Random
using PyCall
using Functors: fmap
using DLPack
using CUDA

if !ispysetup[]
    return
end

fexp = pyimport("functorch.experimental")

if CUDA.functional()
    device = torch.device("cuda:0")
else
    device = torch.device("cpu")
end

model = torch.hub.load("pytorch/vision", "resnet18", pretrained=true)
model_gn = fexp.replace_all_batch_norm_modules_(model).to(device=device)
modelwrap = TorchModuleWrapper(model_gn)
if CUDA.functional()
    modelwrap = fmap(CUDA.cu, modelwrap)
end
x = randn(Float32, reverse((1, 3, 224, 224)))
if CUDA.functional()
    x = CUDA.cu(x)
end

grad,  = Zygote.gradient(m->sum(m(x)), modelwrap)
@test length(grad.params) == length(modelwrap.params)
params = map(x ->  DLPack.share(x, PyObject, pyfrom_dlpack).to(device = device, dtype = modelwrap.dtype).requires_grad_(true), modelwrap.params)
torch_out = modelwrap.torch_stateless_module(params, modelwrap.buffers, map(z-> DLPack.share(z, PyObject, pyfrom_dlpack).to(dtype=modelwrap.dtype, device=device), [x])...).sum()
torchgrad = map(x-> x.cpu().numpy(), torch.autograd.grad(torch_out, params))
@test length(torchgrad) == length(grad.params)
for i in 1:length(grad.params)
    @test isapprox(sum(torchgrad[i]), sum(grad.params[i]), atol=1e-3, rtol=1e-3)
end


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")
include("circular_bc_transformation_elasticity_solver.jl")

function compute_L2_stress_error(polyorder, nelmts, width)
    K1, K2 = 247.0, 192.0    # Pa
    mu1, mu2 = 126.0, 87.0   # Pa
    lambda1 = lame_lambda(K1, mu1)
    lambda2 = lame_lambda(K2, mu2)
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)

    theta0 = -0.067
    stiffness = CutCellDG.HookeStiffness(lambda1, mu1, lambda2, mu2)
    transfstress =
        CutCellDG.plane_strain_transformation_stress(lambda1, mu1, theta0)

    meshwidth = [width, width]
    numqp = required_quadrature_order(polyorder) + 2
    penaltyfactor = 1e2

    dx = width / nelmts
    penalty = penaltyfactor / dx * 0.5 * (lambda1 + mu1 + lambda2 + mu2)

    elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
    levelsetbasis = HermiteTensorProductBasis(2)
    quad = tensor_product_quadrature(2, 4)
    dim, nf = size(interpolation_points(levelsetbasis))
    refpoints = interpolation_points(elasticitybasis)

    interfacecenter = [0.5, 0.5]
    interfaceradius = minimum(meshwidth) / 3.0
    outerradius = 2.0
    analyticalsolution = AnalyticalSolution(
        interfaceradius,
        outerradius,
        interfacecenter,
        lambda1,
        mu1,
        lambda2,
        mu2,
        theta0,
    )

    mesh, cellquads, facequads, interfacequads = construct_mesh_and_quadratures(
        meshwidth,
        nelmts,
        elasticitybasis,
        interfacecenter,
        interfaceradius,
        numqp,
    )
    nodaldisplacement = nodal_displacement(
        mesh,
        elasticitybasis,
        cellquads,
        facequads,
        interfacequads,
        stiffness,
        theta0,
        analyticalsolution,
        penalty,
    )

    stresserror = stress_L2_error(
        nodaldisplacement,
        elasticitybasis,
        cellquads,
        stiffness,
        transfstress,
        theta0,
        mesh,
        x -> core_stress(analyticalsolution),
        x -> shell_stress(analyticalsolution, x),
    )

    den = integral_norm_on_mesh(
        x -> exact_stress(analyticalsolution, x),
        cellquads,
        mesh,
        4,
    )

    normalizedstresserr = stresserror ./ den

    return normalizedstresserr
end

function test_quadratic_elements()
    width = 1.0
    polyorder = 2
    powers = [3, 4, 5]
    nelmts = [2^p + 1 for p in powers]

    err = [compute_L2_stress_error(polyorder, ne, width) for ne in nelmts]

    stresserror = Array(transpose(hcat(err...)))
    dx = width ./ nelmts
    stressrates = mapslices(x -> convergence_rate(dx, x), stresserror, dims = 1)

    @test all(stressrates .> 1.8)
end

function test_cubic_elements()
    width = 1.0
    polyorder = 3
    powers = [3, 4, 5]
    nelmts = [2^p + 1 for p in powers]

    err = [compute_L2_stress_error(polyorder, ne, width) for ne in nelmts]

    dx = width ./ nelmts
    stresserror = Array(transpose(hcat(err...)))

    stressrates = mapslices(x -> convergence_rate(dx, x), stresserror, dims = 1)
    @test all(stressrates .> 2.7)
end

test_quadratic_elements()
test_cubic_elements()


text:
export is_braket
export Ket, Bra

"""
    is_braket( ψ :: Vector; atol=ATOL :: Float64) :: Bool

    is_braket(
        ψ :: Adjoint{T, Vector{T}};
        atol=ATOL :: Float64
    ) where T <: Number :: Bool

Returns `true` if vector `ψ` is a valid bra (or ket):

* `ψ` is a real or complex-valued vector.
* `ψ` is normalized with respect to the bra-ket inner prodcut (`ψ' * ψ == 1`).
"""
function is_braket(ψ :: Vector{<:Number}; atol=ATOL :: Float64)
    isapprox( norm(ψ,2), 1, atol=atol )
end
function is_braket(ψ :: Adjoint{T, Vector{T}}; atol=ATOL :: Float64) where T <: Number
    isapprox( norm(ψ,2), 1, atol=atol )
end

"""
    _braket_error(ψ)

Throws the `DomainError`, `"ψ must be a vector that satisfies \`norm(ψ,2) != 1\`."`.
"""
_braket_error(ψ) = throw(DomainError(ψ, "ψ must be a vector that satisfies `norm(ψ,2) != 1`."))

"""
    Ket(ψ :: Vector{T}; atol=ATOL :: Float64) where T <: Number

A normalized column vector on a complex-valued Hilbert space. If called with
an adjoint vector `ψ'`, the adjoint will be taken on construction.

```julia
Ket(ψ :: Adjoint{Vector{T}}; atol=ATOL :: Float64) where T <: Number
```

Similarly a `Bra` can be converted into a `Ket` via the adjoint.

```julia
Ket(ψ :: Bra{T}; atol=ATOL :: Float64) where T <: Number
```

`Ket(ψ)`, throws a `DomainError` if `ψ` is not normalized.
"""
struct Ket{T<:Number} <: AbstractVector{Number}
    ψ :: Vector{T}
    atol :: Float64
    Ket(
        ψ :: Vector{T}; atol=ATOL :: Float64
    ) where T <: Number = is_braket(ψ, atol=atol) ? new{T}(ψ, atol) : _braket_error(ψ)
    Ket(
        ψ :: Adjoint{T,Vector{T}}; atol=ATOL :: Float64
    ) where T <: Number = is_braket(ψ, atol=atol) ? new{T}(
        ψ', atol
    ) : _braket_error(ψ)
end
Base.size(ket::Ket) = size(ket.ψ)
Base.getindex(ket::Ket, I::Int) = getindex(ket.ψ, I...)
Base.setindex!(ket::Ket, val, I::Int) = (ket.ψ[I...] = val)
is_braket(::Ket) = true

"""
    Bra(ψ :: Adjoint{T,Vector{T}}; atol=ATOL :: Float64) where T <: Number

A row vector on a complex valued HIlbet space. Bras are dual to Kets via the adjoint
If called with an `Vector{<:Number}` type, the adjoint will automatically be taken `ψ'`.

```julia
Bra(ψ :: Vector{<:Number}; atol=ATOL :: Float64)
```

Similarly a `Ket` can be converted into a `Bra` via the adjoint.

```julia
Bra(ψ :: Ket; atol=ATOL :: Float64)
```

`Bra(ψ')`, throws a `DomainError` if `ψ'` is not normalized.
"""
struct Bra{T<:Number} <: AbstractMatrix{Number}
    ψ :: Adjoint{T,Vector{T}}
    atol :: Float64
    Bra(
        ψ :: Vector{T}; atol=ATOL :: Float64
    ) where T <: Number = is_braket(ψ, atol=atol) ? new{T}(ψ', atol) : _braket_error(ψ)
    Bra(
        ψ :: Adjoint{T,Vector{T}}; atol=ATOL :: Float64
    ) where T <: Number = is_braket(ψ, atol=atol) ? new{T}(ψ, atol) : _braket_error(ψ)
end
Base.size(bra::Bra) = size(bra.ψ)
Base.getindex(bra::Bra, I::Vararg{Int,2}) = getindex(bra.ψ, I...)
Base.setindex!(bra::Bra, val, I::Vararg{Int,2}) = (bra.ψ[I...] = val)
is_braket(::Bra) = true

# Bra-Ket conversions
Ket(bra :: Bra) = Ket(bra.ψ', atol=bra.atol)
Bra(ket :: Ket) = Bra(ket.ψ', atol=ket.atol)

"""
Outer product ``|\\psi \\rangle\\langle \\psi|``:

```julia
*(ket :: Ket, bra :: Bra) :: Matrix
```

Inner product ``\\langle \\psi | \\psi \\rangle``:

```julia
*(bra :: Bra, ket :: Ket) :: Number
```
"""
*(ket :: Ket, bra :: Bra) :: Matrix = ket.ψ * bra.ψ
*(bra :: Bra, ket :: Ket) :: Number = (bra.ψ * ket.ψ)[1]

"""
The adjoint is the complex conjugate transpose. This operation converts a `Ket` to
a `Bra`, ``\\langle \\psi | = |\\psi \\rangle^{\\dagger}``:

```julia
adjoint(ket :: Ket) :: Bra
```

Or, convert a `Bra` to a `Ket`, ``|\\psi\\rangle = \\langle \\psi |^{\\dagger}``:

```julia
adjoint(bra :: Bra) :: Ket
```
"""
adjoint(ket :: Ket) = Bra(ket)
adjoint(bra :: Bra) = Ket(bra)

"""
Performs the kronecker product of a set of `Ket`s (`Bra`s) to produce a new `Ket` (`Bra`).

```julia
kron( kets :: Vararg{Ket}; atol=ATOL ) :: Ket
kron( bras :: Vararg{Bra}; atol=ATOL ) :: Bra
```
"""
kron(kets :: Vararg{Ket}; atol=ATOL) = Ket(kron(map(ket -> ket.ψ, kets)...), atol=atol)
kron(bras :: Vararg{Bra}; atol=ATOL) = Bra(kron(map(bra -> bra.ψ', bras)...), atol=atol)


text:
function moonbar()
    try
        print("\x1b[?25l") # hide cursor
        while true
            for 🌝 ∈ '🌑':'🌘'
                sleep(0.1)
                print(🌝)
                print('\r')
            end
        end
    finally
        print("\x1b[?25h") # unhide cursor
    end
end

moonbar()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using PkgDeps
using Test
using UUIDs

const DEPOT = joinpath(@__DIR__, "resources")
const GENERAL_REGISTRY = only(reachable_registries("General"; depots=DEPOT))
const FOOBAR_REGISTRY = only(reachable_registries("Foobar"; depots=DEPOT))

@testset "internal functions" begin
    @testset "_get_pkg_name" begin
        @testset "uuid to name" begin
            expected = "Case1"
            pkg_name = PkgDeps._get_pkg_name(UUID("00000000-1111-2222-3333-444444444444"); depots=DEPOT)

            @test expected == pkg_name
        end

        @testset "exception" begin
            @test_throws NoUUIDMatch PkgDeps._get_pkg_name(UUID("00000000-0000-0000-0000-000000000000"); registries=FOOBAR_REGISTRY)
        end
    end

    @testset "_get_pkg_uuid" begin
        @testset "name to uuid" begin
            expected = UUID("00000000-1111-2222-3333-444444444444")

            pkg_uuid = PkgDeps._get_pkg_uuid("Case1", "Foobar"; depots=DEPOT)
            @test expected == pkg_uuid

            pkg_uuid = PkgDeps._get_pkg_uuid("Case1", FOOBAR_REGISTRY)
            @test expected == pkg_uuid
        end

        @testset "exception" begin
            @test_throws PackageNotInRegistry PkgDeps._get_pkg_uuid("PkgDepsFakePackage", "General")
            @test_throws PackageNotInRegistry PkgDeps._get_pkg_uuid("FakePackage", FOOBAR_REGISTRY)
        end
    end

    @testset "_get_latest_version" begin
        expected = v"0.2.0"
        path = joinpath("resources", "registries", "General", "Case4")
        result = PkgDeps._get_latest_version(path)

        @test expected == result
    end

    @testset "_find_alternative_packages" begin
        MAX = 9
        pkg_to_compare = "package_name"
        packages = ["$(pkg_to_compare)_$(i)" for i in 1:MAX]

        result = PkgDeps._find_alternative_packages(pkg_to_compare, packages)
        @test length(result) == MAX
    end
end

@testset "reachable_registries" begin
    @testset "specfic registry -- $(typeof(v))" for v in ("Foobar", ["Foobar"])
        registry = only(reachable_registries("Foobar"; depots=DEPOT))

        @test registry.name == "Foobar"
    end

    @testset "all registries" begin
        registries = reachable_registries(; depots=DEPOT)

        @test length(registries) == 2
    end
end

@testset "users" begin
    all_registries = reachable_registries(; depots=DEPOT)

    @testset "specific registry - registered in another registry" begin
        dependents = users("DownDep", FOOBAR_REGISTRY; registries=[GENERAL_REGISTRY], depots=DEPOT)

        @test length(dependents) == 1
        [@test case in dependents for case in ["Case3"]]
    end

    @testset "all registries" begin
        dependents = users("DownDep", FOOBAR_REGISTRY; registries=all_registries, depots=DEPOT)

        @test length(dependents) == 3
        @test !("Case4" in dependents)
        [@test case in dependents for case in ["Case1", "Case2", "Case3"]]
    end
end


text:
################################################
#       TEST DE FORWARD PROPAGATION
################################################
using Test

@testset "forward_propagation correct types" begin
    ## Comprobación de tipos y dimensión
    entry_dimesion = 2
    number_of_hide_units = 3
    output_dimension = 2
    OLNN = RandomWeightsNN(
        entry_dimesion,
        number_of_hide_units,
        output_dimension
        )
    # El resultado debe de ser un vector
    @test typeof(forward_propagation(OLNN,ReLU, [1,2.0])) == Vector{Float64}
    # La evaluación debe de tener las mismas dimensiones que la salida de la red neuronal
    @test length(forward_propagation(OLNN,ReLU, [0,1.0])) == output_dimension
end
@testset "forward_propagation matrix order and " begin
    ## Comprobación de evaluación correcta
    # Debiera de ser la red neurona identidad
    S = [0, 0]
    A = [1 0; 0 1]
    B = [1 0; 0 1]
    h = FromMatrixNN(S,A,B)

    vectores = [
        [1,2], [0,0],[-1,4]
    ]
    for v in vectores
        @test forward_propagation(h, x->x,v ) == v
    end
    # Debiera de ser la red neuronal que multiplica el primer índice por dos y el resto por 3
    S = [0, 0]
    A = [1 0; 0 1]
    B = [2 0; 0 3]
    h = FromMatrixNN(S,A,B)

    vectores = [
        [1,2], [0,0],[-1,4]
    ]
    for v in vectores
        @test forward_propagation(h, x->x,v ) == [2*v[1], 3*v[2]]
    end
    # Funcionamiento correcto para translaciones
    # Debiera de ser la red neuronal que suma el vector (1 2)
    S = [1, 2]
    A = [1 0; 0 1]
    B = [1 0; 0 1]
    h = FromMatrixNN(S,A,B)

    vectores = [
        [1,2], [0,0],[-1,4]
    ]
    for v in vectores
        @test forward_propagation(h, x->x,v ) == [v[1]+1, v[2]+2]
    end
    S = [2, 5]
    A = [2 3; 7 8]
    B = [4 7; 10 -9]
    v = [2, -1]
    # Calculamos manualmente cuál debiera de ser el resultado
    c = A*v
    c = c + S
    c = B*c

    h = FromMatrixNN(S,A,B)
    @test forward_propagation(h, x->x,v) == c

end
@testset "forward_propagation activation function" begin
    # Comprobamos que admite una función de activación cualquiera
    S = [0, 0]
    A = [1 0; 0 1]
    B = [1 0; 0 1]
    h = FromMatrixNN(S,A,B)

    vectores = [
        [1,2], [0,-3]
    ]
    soluciones_reLU = [
        [1,2],[0,0]
    ]
    for (v,test) in zip(vectores,soluciones_reLU)
        @test forward_propagation(h, ReLU,v ) == test
    end
    # Comprobamos que aplica correctamente los coeficientes
    S = [0, 0]
    A = [1 0; 0 1]
    B = [1 0; 0 1]
    h = FromMatrixNN(S,A,B)

    vectores = [
        [-1,2], [0,-3]
    ]
    soluciones= [
        [1,4],[0,9]
    ]
    for (v,test) in zip(vectores,soluciones)
        @test forward_propagation(h, x-> x^2,v ) == test
    end

end


text:
module Order

## notions of element ordering ##

export # not exported by Base
    Ordering, Forward, Reverse, Lexicographic,
    By, Lt, Perm,
    ReverseOrdering, ForwardOrdering, LexicographicOrdering,
    DirectOrdering,
    lt, uint_mapping, ord, ordtype

abstract Ordering

immutable ForwardOrdering <: Ordering end
immutable ReverseOrdering{Fwd<:Ordering} <: Ordering
    fwd::Fwd
end

ReverseOrdering(rev::ReverseOrdering) = rev.fwd
ReverseOrdering{Fwd}(fwd::Fwd) = ReverseOrdering{Fwd}(fwd)

typealias DirectOrdering Union(ForwardOrdering,ReverseOrdering{ForwardOrdering})

const Forward = ForwardOrdering()
const Reverse = ReverseOrdering(Forward)

immutable LexicographicOrdering <: Ordering end
const Lexicographic = LexicographicOrdering()

immutable By <: Ordering
    by::Function
end

immutable Lt <: Ordering
    lt::Function
end

immutable Perm{O<:Ordering,V<:AbstractVector} <: Ordering
    order::O
    data::V
end
Perm{O<:Ordering,V<:AbstractVector}(o::O,v::V) = Perm{O,V}(o,v)

lt(o::ForwardOrdering,       a, b) = isless(a,b)
lt(o::ReverseOrdering,       a, b) = lt(o.fwd,b,a)
lt(o::By,                    a, b) = isless(o.by(a),o.by(b))
lt(o::Lt,                    a, b) = o.lt(a,b)
lt(o::LexicographicOrdering, a, b) = lexcmp(a,b) < 0

function lt(p::Perm, a::Int, b::Int)
    lt(p.order, p.data[a], p.data[b]) ? true :
    lt(p.order, p.data[b], p.data[a]) ? false : a < b
end
function lt(p::Perm{LexicographicOrdering}, a::Int, b::Int)
    c = lexcmp(p.data[a], p.data[b])
    c != 0 ? c < 0 : a < b
end

# Map a bits-type to an unsigned int, maintaining sort order
uint_mapping(::ForwardOrdering, x::Unsigned) = x
uint_mapping(::ForwardOrdering, x::Int8)     = uint8  (x $ typemin(Int8))
uint_mapping(::ForwardOrdering, x::Int16)    = uint16 (x $ typemin(Int16))
uint_mapping(::ForwardOrdering, x::Int32)    = uint32 (x $ typemin(Int32))
uint_mapping(::ForwardOrdering, x::Int64)    = uint64 (x $ typemin(Int64))
uint_mapping(::ForwardOrdering, x::Int128)   = uint128(x $ typemin(Int128))
uint_mapping(::ForwardOrdering, x::Float32)  = (y = reinterpret(Int32, x); uint32(y < 0 ? ~y : (y $ typemin(Int32))))
uint_mapping(::ForwardOrdering, x::Float64)  = (y = reinterpret(Int64, x); uint64(y < 0 ? ~y : (y $ typemin(Int64))))

uint_mapping{Fwd}(::ReverseOrdering{Fwd}, x) = ~uint_mapping(Fwd, x)
#uint_mapping{T<:Real}(::ReverseOrdering{ForwardOrdering}, x::T) = ~uint_mapping(Forward, x)  ## Manually inlined
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Unsigned) = ~x
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Int8)     = ~uint8  (x $ typemin(Int8))
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Int16)    = ~uint16 (x $ typemin(Int16))
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Int32)    = ~uint32 (x $ typemin(Int32))
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Int64)    = ~uint64 (x $ typemin(Int64))
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Int128)   = ~uint128(x $ typemin(Int128))
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Float32)  = (y = reinterpret(Int32, x); uint32(y < 0 ? y : ~(y $ typemin(Int32))))
uint_mapping(::ReverseOrdering{ForwardOrdering}, x::Float64)  = (y = reinterpret(Int64, x); uint64(y < 0 ? y : ~(y $ typemin(Int64))))

uint_mapping(o::By,   x     ) = uint_mapping(Forward, o.by(x))
uint_mapping(o::Perm, i::Int) = uint_mapping(o.order, o.data[i])
uint_mapping(o::Lt,   x     ) = error("uint_mapping does not work with general Lt Orderings")

ordtype(o::ReverseOrdering, vs::AbstractArray) = ordtype(o.fwd, vs)
ordtype(o::Perm,            vs::AbstractArray) = ordtype(o.order, o.data)
# TODO: here, we really want the return type of o.by, without calling it
ordtype(o::By,              vs::AbstractArray) = try typeof(o.by(vs[1])) catch Any end
ordtype(o::Ordering,        vs::AbstractArray) = eltype(vs)

function ord(lt::Function, by::Function, rev::Bool, order::Ordering=Forward)
    order == Forward || order == Lexicographic ||
        Base.warn_once("the `order` keyword is deprecated, use `lt`, `by` and `rev` instead")
    o = (lt===isless) & (by===identity) ? order  :
        (lt===isless) & (by!==identity) ? By(by) :
        (lt!==isless) & (by===identity) ? Lt(lt) :
                                          Lt((x,y)->lt(by(x),by(y)))
    rev ? ReverseOrdering(o) : o
end

end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
# no kwargs
precompile(calendar, ())
precompile(calendar, (Date, ))

precompile(calendar, (Vector{Date}, ))
precompile(calendar, (StepRange{Date, Month}, ))

precompile(calendar, (Int, ))
precompile(calendar, (Int, Int))

precompile(calendar, (Vector{Int64}, Int))
precompile(calendar, (UnitRange{Int64}, Int))
precompile(calendar, (StepRange{Int64}, Int))

precompile(calendar, (Int, Vector{Int64}))
precompile(calendar, (Int, UnitRange{Int64}))
precompile(calendar, (Int, StepRange{Int64}))

precompile(calendar, (Vector{Int64}, Vector{Int64}))
precompile(calendar, (Vector{Int64}, UnitRange{Int64}))
precompile(calendar, (Vector{Int64}, StepRange{Int64}))

precompile(calendar, (UnitRange{Int64}, Vector{Int64}))
precompile(calendar, (UnitRange{Int64}, UnitRange{Int64}))
precompile(calendar, (UnitRange{Int64}, StepRange{Int64}))

precompile(calendar, (StepRange{Int64}, Vector{Int64}))
precompile(calendar, (StepRange{Int64}, UnitRange{Int64}))
precompile(calendar, (StepRange{Int64}, StepRange{Int64}))

# with kwargs
precompile(calendar, (TextFormat,))
precompile(calendar, (Date, TextFormat))

precompile(calendar, (Vector{Date}, TextFormat))
precompile(calendar, (StepRange{Date, Month}, TextFormat))

precompile(calendar, (Int, TextFormat))
precompile(calendar, (Int, Int, TextFormat))

precompile(calendar, (Vector{Int64}, Int, TextFormat))
precompile(calendar, (UnitRange{Int64}, Int, TextFormat))
precompile(calendar, (StepRange{Int64}, Int, TextFormat))

precompile(calendar, (Int, Vector{Int64}, TextFormat))
precompile(calendar, (Int, UnitRange{Int64}, TextFormat))
precompile(calendar, (Int, StepRange{Int64}, TextFormat))

precompile(calendar, (Vector{Int64}, Vector{Int64}, TextFormat))
precompile(calendar, (Vector{Int64}, UnitRange{Int64}, TextFormat))
precompile(calendar, (Vector{Int64}, StepRange{Int64}, TextFormat))

precompile(calendar, (UnitRange{Int64}, Vector{Int64}, TextFormat))
precompile(calendar, (UnitRange{Int64}, UnitRange{Int64}, TextFormat))
precompile(calendar, (UnitRange{Int64}, StepRange{Int64}, TextFormat))

precompile(calendar, (StepRange{Int64}, Vector{Int64}, TextFormat))
precompile(calendar, (StepRange{Int64}, UnitRange{Int64}, TextFormat))
precompile(calendar, (StepRange{Int64}, StepRange{Int64}, TextFormat))


text:
export logsumexp, rDirichlet, embed, condNorm,
ldnorm, lmvbeta, vech, xpnd_tri;

"""
    logsumexp(x[, usemax])

Computes `log(sum(exp(x)))` in a stable manner.

### Example
```julia
  x = rand(5)
  logsumexp(x)
  log(sum(exp.(x)))
```
"""
function logsumexp(x::Array{T}, usemax::Bool=true) where T <: Real
  if usemax
    m = maximum(x)
  else
    m = minimum(x)
  end

  m + log(sum(exp.(x .- m)))
end

"""
    logsumexp(x, region[, usemax])

Computes `log(sum(exp(x)))` in a stable manner along dimensions specified.

### Example
```julia
  x = reshape(collect(1:24)*1.0, (2,3,4))
  logsumexp(x, 2)
```
"""
function logsumexp(x::Array{T}, region, usemax::Bool=true) where T <: Real
  if usemax
    ms = maximum(x, dims=region)
  else
    ms = minimum(x, dims=region)
  end
  bc_xminusms = broadcast(-, x, ms)

  expxx = exp.(bc_xminusms)
  sumexpxx = sum(expxx, dims=region)

  log.(sumexpxx) .+ ms
end

"""
    rDirichlet(α[, logscale])

  Single draw from Dirichlet distribution, option for log scale.

  ### Example
  ```julia
  rDirichlet(ones(5),true)
  ```
"""
function rDirichlet(α::Array{T, 1}, logscale::Bool=false) where T <: Real
  @assert( all(α .> 0.0) )

  k = length(α)
  xx = Vector{T}(undef, k) # allows changes to elements
  s = 0.0

  if logscale
    for i in 1:k
      xx[i] = log(rand(Gamma(α[i], 1.0)))
    end
    s = logsumexp(xx)
    out = xx .- s

  else
    for i in 1:k
      xx[i] = rand(Gamma(α[i], 1.0))
      s += xx[i]
    end
    out = xx / s

  end

out
end

"""
    embed(y::Union{Vector{Int}, Vector{T}}, nlags::Int)

### Example
```julia
embed(collect(1:5), 2)
```
"""
function embed(y::Union{Vector{Int}, Vector{Float64}}, nlags::Int)
    TT = length(y)
    n = TT - nlags
    emdim = nlags + 1
    out = Array{Float64, 2}(undef, n, emdim)
    for i in 1:n
        tt = i + nlags
        out[i,:] = copy(y[range(tt, step=-1, length=emdim)])
    end
    out
end

"""
    condNorm(μx, μy, Σxx, Σxy, Σyy, x)

Returns the conditional distribution of a the ``y`` subvector
  of ``(x,y) ∼ MultivariateNormal``.

### Example
```julia
  xx = randn(10, 5)
  Sig = xx'xx
  mu = collect(1:5)*1.0
  x_indx = 1:3
  y_indx = 4:5
  μx = mu[x_indx]
  μy = mu[y_indx]
  Σxx = Sig[x_indx, x_indx]
  Σyy = Sig[y_indx, y_indx]
  Σxy = Sig[x_indx, y_indx]
  x = [0.3, 2.1, 3.4]
  @time condNorm(μx, μy, Σxx, Σxy, Σyy, x)
```
"""
function condNorm(μx::Union{T, Vector{T}},
  μy::Union{T, Vector{T}},
  Σxx::Union{T, Array{T, 2}},
  Σxy::Union{T, Vector{T}, Array{T, 2}},
  Σyy::Union{T, Array{T, 2}}, x::Union{T, Vector{T}}) where T <: Real

  isposdef(Σxx) || throw(ArgumentError("Σxx must be positive definite."))
  isposdef(Σyy) || throw(ArgumentError("Σyy must be positive definite."))

  # Cholesky method (about 2 times faster than direct method)
  U = (cholesky(Σxx)).U
  Ut = transpose(U)

  # C = U' \ Σxy
  # C = At_ldiv_B(U, Σxy)
  C = Ut \ Σxy
  # D = U' \ (x - μx)
  # D = At_ldiv_B(U, (x - μx))
  D = Ut \ (x - μx)

  μ = μy + C'D
  Σ = Σyy - C'C

  # Direct method
  # A = Σxx \ (x - μx)
  # B = Σxx \ Σxy
  #
  # μ = μy + Σxy'A
  # Σ = Σyy - Σxy'B
  # Σ = 0.5 * (Σ' + Σ)

  # because univariate normal parameterized with stdev
  if length(μ) == 1
    Σ = sqrt(Σ)
  end

  out = MultivariateNormal(μ, Σ)
end

## triangular matrix from vector
# https://groups.google.com/forum/#!topic/julia-users/UARlZBCNlng
# [ i<=j ? v[j*(j-1)>>>1+i] : 0 for i=1:n, j=1:n ]

"""
    ldnorm(x, μ, σ2)
### Example
```julia
ldnorm(0.5, 0.0, 1.0)
```
"""
function ldnorm(x::T, μ::T, σ2::T) where T <: Real
  -0.5*log(2*π*σ2) - 0.5 * (x - μ)^2 / σ2
end

"""
    lmvbeta(x)

Computes the natural log of ``∏(Γ(x)) / Γ(sum(x))``.

### Example
```julia
  x = rand(5)
  lmvbeta(x)
```
"""
function lmvbeta(x::Array{T}) where T <: Real
    sum(lgamma.(x)) - lgamma(sum(x))
end

"""
    vech(A)

Extracts the lower triangular elements of a matrix as a vector.

Adapted from https://discourse.julialang.org/t/half-vectorization/7399/4

### Example
```julia
    A = reshape(collect(1:16), 4,4)
    v = vech(A)
    xpnd_tri(v)

    v = vech(A, false)
    xpnd_tri(v)

    vech(A', false)
```
"""
function vech(A::AbstractArray{T, 2},
    diag::Bool=true) where {T}
    m = LinearAlgebra.checksquare(A)
    k = 0
    if diag
        v = Vector{T}(undef, (m*(m+1))>>1)
        for j = 1:m, i = j:m
            @inbounds v[k += 1] = A[i,j]
        end
    else
        v = Vector{T}(undef, (m*(m-1))>>1)
        for j = 1:(m-1), i = (j+1):m
            @inbounds v[k += 1] = A[i,j]
        end
    end
    return v
end

"""
    xpnd(v)

Creates a lower triangular matrix from the elements of a vector.

Adapted from https://groups.google.com/forum/#!topic/julia-users/UARlZBCNlng
and xpnd function from R package MCMCpack.

### Example
```julia
  A = reshape(collect(1:16), 4,4)
  v = vech(A)
  xpnd_tri(v)

  v = vech(A, false)
  xpnd_tri(v, false)
```
"""
function xpnd_tri(v::Vector{T}, diag::Bool=true) where {T}
    k = 0
    n = length(v)
    if diag
        m = Int((sqrt(1 + 8*n) - 1) / 2)
        A0 = Array{T, 2}(undef, m, m)
        for j = 1:m, i = j:m
            k += 1
            @inbounds A0[i,j] = v[k]
        end
        A = LowerTriangular(A0)
    else
        m = Int((sqrt(1 + 8*n) + 1) / 2)
        A0 = Array{T, 2}(undef, m, m)
        for j = 1:(m-1), i = (j+1):m
            k += 1
            @inbounds A0[i,j] = v[k]
        end
        A = UnitLowerTriangular(A0)
    end
    return A
end


text:
for (elty, relty, ext) in ((:Float32, :Float32, :s),
                           (:Float64, :Float64, :d))
    for (matA, matB, sym) in ((:Matrix, :Matrix, "_"),
                              (:DistMatrix, :DistMatrix, "Dist_"),
                              (:SparseMatrix, :Matrix, "Sparse_"),
                              (:DistSparseMatrix, :DistMultiVec, "DistSparse_"))
        @eval begin
            function leastSquares!(A::$matA{$elty}, B::$matB{$elty}, X::$matB{$elty};
                orientation::Integer = NORMAL)
                err = ccall(($(string("ElLeastSquares", sym, ext)), libEl), Cuint,
                    (Cuint, Ptr{Void}, Ptr{Void}, Ptr{Void}),
                    orientation, A.obj, B.obj, X.obj)
                err == 0 || throw(ElError(err))
                return X
            end
        end
    end
end

function leastSquares{T}(A::DistMatrix{T}, B::DistMatrix{T};
    orientation::Integer = NORMAL)
    X = DistMatrix(T, MC, MR, Grid(A))
    return leastSquares!(A, B, X, orientation = orientation)
end
function leastSquares{T}(A::DistSparseMatrix{T}, B::DistMultiVec{T};
    orientation::Integer = NORMAL)
    X = DistMultiVec(T, comm(A))
    return leastSquares!(A, B, X, orientation = orientation)
end


text:
using MonteCarloCollisions
using Statistics
using LinearAlgebra
using Test

@testset "MonteCarloCollisions.jl" begin
    nsampl = 1000
    samples = [MonteCarloCollisions.random_direction(1.) for i in 1:nsampl]
    @test max(norm.(samples)...) ≈ min(norm.(samples)...)
    @test mean(samples) < std(samples)/sqrt(nsampl)*5
    # Write your tests here.
end