text:
# # Incompressible Navier-Stokes Equations via DifferentialEquations.jl
#
# ![](https://user-images.githubusercontent.com/9196588/134514213-76d91d34-19ab-47c2-957e-16bb0c8669e1.gif)
#
#
# In this example we focus on a simple but visually appealing problem from
# fluid dynamics, namely vortex shedding. This problem is also known as
# [von-Karman vortex streets](https://en.wikipedia.org/wiki/K%C3%A1rm%C3%A1n_vortex_street). Within this example, we show how to utilize [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl)
# in tandem with Ferrite.jl to solve this space-time problem. To keep things simple we use a naive approach
# to discretize the system.
#
# ## Remarks on DifferentialEquations.jl
#
# Many "time step solvers" of [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl) assume that that the
# problem is provided in mass matrix form. The incompressible Navier-Stokes
# equations as stated above yield a DAE in this form after applying a spatial
# discretization technique - in our case FEM. The mass matrix form of ODEs and DAEs
# is given as:
# ```math
#   M(t) \mathrm{d}_t u = f(u,t)
# ```
# where $M$ is a possibly time-dependent and not necessarily invertible mass matrix,
# $u$ the vector of unknowns and $f$ the right-hand-side (RHS). For us $f$ can be interpreted as
# the spatial discretization of all linear and nonlinear operators depending on $u$ and $t$,
# but not on the time derivative of $u$.
#
# ## Some Theory on the Incompressible Navier-Stokes Equations
#
# ### Problem Description in Strong Form
#
# The incompressible Navier-Stokes equations can be stated as the system
# ```math
#  \begin{aligned}
#    \partial_t v &= \underbrace{\nu \Delta v}_{\text{viscosity}} - \underbrace{(v \cdot \nabla) v}_{\text{advection}} - \underbrace{\nabla p}_{\text{pressure}} \\
#               0 &= \underbrace{\nabla \cdot v}_{\text{incompressibility}}
#  \end{aligned}
# ```
# where $v$ is the unknown velocity field, $p$ the unknown pressure field,
# $\nu$ the dynamic viscosity and $\Delta$ the Laplacian. In the derivation we assumed
# a constant density of 1 for the fluid and negligible coupling between the velocity components.
# Finally we see that the pressure term appears only in combination with the gradient
# operator, so for any solution $p$ the function $p + c$ is also an admissible solution, if
# we do not impose Dirichlet conditions on the pressure. To resolve this we introduce the
# implicit constraint that $ \int_\Omega p = 0 $.
#
# Our setup is derived from [Turek's DFG benchmark](http://www.mathematik.tu-dortmund.de/~featflow/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark1_re20.html).
# We model a channel with size $0.41 \times 2.2$ and a hole of radius $0.05$ centered at $(0.2, 0.2)$.
# The left side has a parabolic inflow profile, which is ramped up over time, modeled as the time dependent
# Dirichlet condition
# ```math
#  v(x,y,t)
#  =
#  \begin{bmatrix}
#      4 v_{in}(t) y (0.41-y)/0.41^2 \\
#      0
#  \end{bmatrix}
# ```
# where $v_{in}(t) = \text{clamp}(t, 0.0, 1.0)$. With a dynamic viscosity of $\nu = 0.001$
# this is enough to induce turbulence behind the cylinder which leads to vortex shedding. The top and bottom of our
# channel have no-slip conditions, i.e. $v = [0,0]^{\textrm{T}}$, while the right boundary has the do-nothing boundary condtion
# $\nu \partial_{\textrm{n}} v - p n = 0$ to model outflow. With these boundary conditions we can choose the zero solution as a
# feasible initial condition.
#
# ### Derivation of Semi-Discrete Weak Form
#
# By multiplying test functions $\varphi$ and $\psi$ from a suitable test function space on the strong form,
# followed by integrating over the domain and applying partial integration to the pressure and viscosity terms
# we can obtain the following weak form
# ```math
#  \begin{aligned}
#    \int_\Omega \partial_t v \cdot \varphi &= - \int_\Omega \nu \nabla v : \nabla \varphi - \int_\Omega (v \cdot \nabla) v \cdot \varphi + \int_\Omega p (\nabla \cdot \varphi) + \int_{\partial \Omega_{N}} \underbrace{(\nu \partial_n v - p n )}_{=0} \cdot \varphi \\
#                                  0 &= \int_\Omega (\nabla \cdot v) \psi
#  \end{aligned}
# ```
# for all possible test functions from the suitable space.
#
# Now we can discretize the problem as usual with the finite element method
# utilizing Taylor-Hood elements (Q2Q1) to yield a stable discretization in
# mass matrix form:
# ```math
#  \underbrace{\begin{bmatrix}
#      M_v & 0 \\
#       0  & 0
#  \end{bmatrix}}_{:=M}
#  \begin{bmatrix}
#      \mathrm{d}_t\hat{v} \\
#      \mathrm{d}_t\hat{p}
#  \end{bmatrix}
#  =
#  \underbrace{\begin{bmatrix}
#       A & B^{\textrm{T}} \\
#       B & 0
#  \end{bmatrix}}_{:=K}
#  \begin{bmatrix}
#      \hat{v} \\
#      \hat{p}
#  \end{bmatrix}
#  +
#  \begin{bmatrix}
#      N(\hat{v}, \hat{v}, \hat{\varphi}) \\
#      0
#  \end{bmatrix}
# ```
# Here $M$ is the singular block mass matrix, $K$ is the discretized Stokes operator and $N$ the nonlinear advection term, which
# is also called trilinear form. $\hat{v}$ and $\hat{p}$ represent the time-dependent vectors of nodal values of the discretizations
# of $v$ and $p$ respectively, while $\hat{\varphi}$ is the choice for the test function in the discretization. The hats are dropped
# in the implementation and only stated for clarity in this section.
#
#
# ## Commented Implementation
#
# Now we solve the problem with Ferrite and [DifferentialEquations.jl](https://github.com/SciML/DifferentialEquations.jl). What follows is a program spliced with comments.
# The full program, without comments, can be found in the next [section](@ref ns_vs_diffeq-plain-program).
#
# First we load Ferrite and some other packages we need
using Ferrite, SparseArrays, BlockArrays, LinearAlgebra, UnPack
# Since we do not need the complete DifferentialEquations suite, we just load the required ODE infrastructure, which can also handle
# DAEs in mass matrix form.
using OrdinaryDiffEq

# We start off by defining our only material parameter.
ν = 1.0/1000.0; #dynamic viscosity

# Next a fine 2D rectangular grid has to be generated. We leave the cell size parametric for flexibility when
# playing around with the code. Note that the mesh is pretty fine, leading to a high memory consumption when
# feeding the equation system to direct solvers.
dim = 2
cell_scale_factor = 2.0
x_cells = round(Int, cell_scale_factor*220)
y_cells = round(Int, cell_scale_factor*41)
# CI chokes if the grid is too fine. :)     #src
x_cells = round(Int, 55/3)                  #hide
y_cells = round(Int, 41/3)                  #hide
grid = generate_grid(Quadrilateral, (x_cells, y_cells), Vec{2}((0.0, 0.0)), Vec{2}((2.2, 0.41)));

# Next we carve a hole $B_{0.05}(0.2,0.2)$ in the mesh by deleting the cells and update the boundary face sets.
# This code will be replaced once a proper mesh interface is avaliable.
cell_indices = filter(ci->norm(mean(map(i->grid.nodes[i].x-[0.2,0.2], Ferrite.vertices(grid.cells[ci]))))>0.05, 1:length(grid.cells))
hole_cell_indices = filter(ci->norm(mean(map(i->grid.nodes[i].x-[0.2,0.2], Ferrite.vertices(grid.cells[ci]))))<=0.05, 1:length(grid.cells));
hole_face_ring = Set{FaceIndex}()
for hci ∈ hole_cell_indices
    push!(hole_face_ring, FaceIndex((hci+1, 4)))
    push!(hole_face_ring, FaceIndex((hci-1, 2)))
    push!(hole_face_ring, FaceIndex((hci-x_cells, 3)))
    push!(hole_face_ring, FaceIndex((hci+x_cells, 1)))
end
grid.facesets["hole"] = Set(filter(x->x.idx[1] ∉ hole_cell_indices, collect(hole_face_ring)));
cell_indices_map = map(ci->norm(mean(map(i->grid.nodes[i].x-[0.2,0.2], Ferrite.vertices(grid.cells[ci]))))>0.05 ? indexin([ci], cell_indices)[1] : 0, 1:length(grid.cells))
grid.cells = grid.cells[cell_indices]
for facesetname in keys(grid.facesets)
    grid.facesets[facesetname] = Set(map(fi -> FaceIndex( cell_indices_map[fi.idx[1]] ,fi.idx[2]), collect(grid.facesets[facesetname])))
end;

# We test against full development of the flow - so regenerate the grid                              #src
grid = generate_grid(Quadrilateral, (x_cells, y_cells), Vec{2}((0.0, 0.0)), Vec{2}((0.55, 0.41)));   #hide

# ### Function Space
# To ensure stability we utilize the Taylor-Hood element pair Q2-Q1.
# We have to utilize the same quadrature rule for the pressure as for the velocity, because in the weak form the
# linear pressure term is tested against a quadratic function.
ip_v = Lagrange{dim, RefCube, 2}()
ip_geom = Lagrange{dim, RefCube, 1}()
qr = QuadratureRule{dim, RefCube}(4)
cellvalues_v = CellVectorValues(qr, ip_v, ip_geom);

ip_p = Lagrange{dim, RefCube, 1}()
cellvalues_p = CellScalarValues(qr, ip_p, ip_geom);

dh = DofHandler(grid)
push!(dh, :v, dim, ip_v)
push!(dh, :p, 1, ip_p)
close!(dh);

# ### Boundary Conditions
# As in the DFG benchmark we apply no-slip conditions to the top, bottom and
# cylinder boundary. The no-slip condition states that the velocity of the
# fluid on this portion of the boundary is fixed to be zero.
ch = ConstraintHandler(dh);

nosplip_face_names = ["top", "bottom", "hole"];
# No hole for the test present                                          #src
nosplip_face_names = ["top", "bottom"]                                  #hide
∂Ω_noslip = union(getfaceset.((grid, ), nosplip_face_names)...);
noslip_bc = Dirichlet(:v, ∂Ω_noslip, (x, t) -> [0,0], [1,2])
add!(ch, noslip_bc);

# The left boundary has a parabolic inflow with peak velocity of 1.0. This
# ensures that for the given geometry the Reynolds number is 100, which
# is already enough to obtain some simple vortex streets. By increasing the
# velocity further we can obtain stronger vortices - which may need additional
# refinement of the grid.
∂Ω_inflow = getfaceset(grid, "left");

vᵢₙ(t) = clamp(t, 0.0, 1.0)*1.0 #inflow velocity
vᵢₙ(t) = clamp(t, 0.0, 1.0)*0.3 #hide
parabolic_inflow_profile((x,y),t) = [4*vᵢₙ(t)*y*(0.41-y)/0.41^2,0]
inflow_bc = Dirichlet(:v, ∂Ω_inflow, parabolic_inflow_profile, [1,2])
add!(ch, inflow_bc);

# The outflow boundary condition has been applied on the right side of the
# cylinder when the weak form has been derived by setting the boundary integral
# to zero. It is also called the do-nothing condition. Other outflow conditions
# are also possible.
∂Ω_free = getfaceset(grid, "right");

close!(ch)
update!(ch, 0.0);

# ### Linear System Assembly
# Next we describe how the block mass matrix and the Stokes matrix are assembled.
#
# For the block mass matrix $M$ we remember that only the first equation had a time derivative
# and that the block mass matrix corresponds to the term arising from discretizing the time
# derivatives. Hence, only the upper left block has non-zero components.
function assemble_mass_matrix(cellvalues_v::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim}, M::SparseMatrixCSC, dh::DofHandler) where {dim}
    ## Allocate a buffer for the local matrix and some helpers, together with the assembler.
    n_basefuncs_v = getnbasefunctions(cellvalues_v)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)
    n_basefuncs = n_basefuncs_v + n_basefuncs_p
    v▄, p▄ = 1, 2
    Mₑ = PseudoBlockArray(zeros(n_basefuncs, n_basefuncs), [n_basefuncs_v, n_basefuncs_p], [n_basefuncs_v, n_basefuncs_p])

    ## It follows the assembly loop as explained in the basic tutorials.
    mass_assembler = start_assemble(M)
    @inbounds for cell in CellIterator(dh)
        fill!(Mₑ, 0)
        Ferrite.reinit!(cellvalues_v, cell)

        for q_point in 1:getnquadpoints(cellvalues_v)
            dΩ = getdetJdV(cellvalues_v, q_point)
            ## Remember that we assemble a vector mass term, hence the dot product.
            for i in 1:n_basefuncs_v
                φᵢ = shape_value(cellvalues_v, q_point, i)
                for j in 1:n_basefuncs_v
                    φⱼ = shape_value(cellvalues_v, q_point, j)
                    Mₑ[BlockIndex((v▄, v▄), (i, j))] += φᵢ ⋅ φⱼ * dΩ
                end
            end
        end
        assemble!(mass_assembler, celldofs(cell), Mₑ)
    end

    return M
end;

# Next we discuss the assembly of the Stokes matrix.
# Remember that we use the same function spaces for trial and test, hence the
# matrix has the following block form
# ```math
#   K = \begin{bmatrix}
#       A & B^{\textrm{T}} \\
#       B & 0
#   \end{bmatrix}
# ```
# which is also called saddle point matrix. These problems are known to have
# a non-trivial kernel, which is a reflection of the strong form as discussed
# in the theory portion if this example.
function assemble_stokes_matrix(cellvalues_v::CellVectorValues{dim}, cellvalues_p::CellScalarValues{dim}, ν, K::SparseMatrixCSC, dh::DofHandler) where {dim}
    ## Again, some buffers and helpers
    n_basefuncs_v = getnbasefunctions(cellvalues_v)
    n_basefuncs_p = getnbasefunctions(cellvalues_p)
    n_basefuncs = n_basefuncs_v + n_basefuncs_p
    v▄, p▄ = 1, 2
    Kₑ = PseudoBlockArray(zeros(n_basefuncs, n_basefuncs), [n_basefuncs_v, n_basefuncs_p], [n_basefuncs_v, n_basefuncs_p])

    ## Assembly loop
    stiffness_assembler = start_assemble(K)
    @inbounds for cell in CellIterator(dh)
        ## Don't forget to initialize everything
        fill!(Kₑ, 0)

        Ferrite.reinit!(cellvalues_v, cell)
        Ferrite.reinit!(cellvalues_p, cell)

        for q_point in 1:getnquadpoints(cellvalues_v)
            dΩ = getdetJdV(cellvalues_v, q_point)
            # Assemble local viscosity block of $A$
            #+
            for i in 1:n_basefuncs_v
                ∇φᵢ = shape_gradient(cellvalues_v, q_point, i)
                for j in 1:n_basefuncs_v
                    ∇φⱼ = shape_gradient(cellvalues_v, q_point, j)
                    Kₑ[BlockIndex((v▄, v▄), (i, j))] -= ν * ∇φᵢ ⊡ ∇φⱼ * dΩ
                end
            end
            # Assemble local pressure and incompressibility blocks of $B^{\textrm{T}}$ and $B$.
            #+
            for j in 1:n_basefuncs_p
                ψ = shape_value(cellvalues_p, q_point, j)
                for i in 1:n_basefuncs_v
                    divφ = shape_divergence(cellvalues_v, q_point, i)
                    Kₑ[BlockIndex((v▄, p▄), (i, j))] += (divφ * ψ) * dΩ
                    Kₑ[BlockIndex((p▄, v▄), (j, i))] += (ψ * divφ) * dΩ
                end
            end
        end

        ## Assemble `Kₑ` into the Stokes matrix `K`.
        assemble!(stiffness_assembler, celldofs(cell), Kₑ)
    end
    return K
end;

# ### Solution of the semi-discretized system via DifferentialEquations.jl
# First we assemble the linear portions for efficiency. These matrices are
# assumed to be constant over time.
T = 10.0
Δt₀ = 0.01
Δt_save = 0.1

M = create_sparsity_pattern(dh);
M = assemble_mass_matrix(cellvalues_v, cellvalues_p, M, dh);

K = create_sparsity_pattern(dh);
K = assemble_stokes_matrix(cellvalues_v, cellvalues_p, ν, K, dh);

# These are our initial conditions. We start from the zero solution, because it
# is trivially admissible if the Dirichlet conditions are zero everywhere on the
# Dirichlet boundary for $t=0$. Note that the time stepper is also doing fine if the
# Dirichlet condition is non-zero and not too pathological.
u₀ = zeros(ndofs(dh))
apply!(u₀, ch);

# DifferentialEquations assumes dense matrices by default, which is not
# feasible for semi-discretization of finite element models. We communicate
# that a sparse matrix with specified pattern should be utilized through the
# `jac_prototyp` argument. It is simple to see that the Jacobian and the
# stiffness matrix share the same sparsity pattern, since they share the
# same relation between trial and test functions.
jac_sparsity = sparse(K);

# To apply the nonlinear portion of the Navier-Stokes problem we simply hand
# over the dof handler and cell values to the right-hand-side (RHS) as a parameter.
# Further the pre-assembled linear part (which is time independent) is
# passed to save some runtime. To apply the time-dependent Dirichlet BCs, we
# also hand over the constraint handler.
# The basic idea to apply the Dirichlet BCs consistently is that we copy the
# current solution `u`, apply the Dirichlet BCs on the copy, evaluate the
# discretized RHS of the Navier-Stokes equations with this vector
# and finally set the RHS to zero on every constraint. This way we obtain a
# correct solution for all dofs which are not Dirichlet constrained. These
# dofs are then corrected in a post-processing step, when evaluating the
# solution vector at specific time points.
# It should be finally noted that this **trick does not work** out of the box
# **for constraining algebraic portion** of the DAE, i.e. if we would like to
# put a Dirichlet BC on pressure dofs. As a workaround we have to set $f_{\textrm{i}} = 1$
# instead of $f_{\textrm{i}} = 0$, because otherwise the equation system gets singular.
# This is obvious when we remember that our mass matrix is zero for these
# dofs, such that we obtain the equation $0 \cdot \mathrm{d}_t p_{\textrm{i}} = 1 \cdot p_{\textrm{i}}$, which
# now has a unique solution.
struct RHSparams
    K::SparseMatrixCSC
    ch::ConstraintHandler
    dh::DofHandler
    cellvalues_v::CellVectorValues
end
p = RHSparams(K, ch, dh, cellvalues_v)

function navierstokes!(du,u_uc,p,t)
    # Unpack the struct to save some allocations.
    #+
    @unpack K,ch,dh,cellvalues_v = p

    # We start by applying the time-dependent Dirichlet BCs. Note that we are
    # not allowed to mutate `u_uc`! We also can not pre-allocate this variable
    # if we want to use AD to derive the Jacobian matrix, which appears in the
    # utilized implicit Euler. If we hand over the Jacobian analytically to
    # the solver, or when utilizing a method which does not require building the
    # Jacobian, then we could also hand over a buffer for `u` in our RHSparams
    # structure to save the allocations made here.
    #+
    u = copy(u_uc)
    update!(ch, t)
    apply!(u, ch)

    # Now we apply the rhs of the Navier-Stokes equations
    #+
    ## Linear contribution (Stokes operator)
    mul!(du, K, u) # du .= K * u

    ## nonlinear contribution
    n_basefuncs = getnbasefunctions(cellvalues_v)
    for cell in CellIterator(dh)
        Ferrite.reinit!(cellvalues_v, cell)
        all_celldofs = celldofs(cell)
        v_celldofs = all_celldofs[dof_range(dh, :v)]
        v_cell = u[v_celldofs]
        for q_point in 1:getnquadpoints(cellvalues_v)
            dΩ = getdetJdV(cellvalues_v, q_point)
            ∇v = function_gradient(cellvalues_v, q_point, v_cell)
            v = function_value(cellvalues_v, q_point, v_cell)
            for j in 1:n_basefuncs
                φⱼ = shape_value(cellvalues_v, q_point, j)
                # Note that in Tensors.jl the definition $\textrm{grad} v = \nabla v$ holds.
                # With this information it can be quickly shown in index notation that
                # ```math
                # [(v \cdot \nabla) v]_{\textrm{i}} = v_{\textrm{j}} (\partial_{\textrm{j}} v_{\textrm{i}}) = [v (\nabla v)^{\textrm{T}}]_{\textrm{i}}
                # ```
                # where we should pay attentation to the transpose of the gradient.
                #+
                du[v_celldofs[j]] -= v ⋅ ∇v' ⋅ φⱼ * dΩ
            end
        end
    end

    # For now we have to ingore the evolution of the Dirichlet BCs.
    # The DBC dofs in the solution vector will be corrected in a post-processing step.
    #+
    apply_zero!(du, ch)
end;
# Finally, together with our pre-assembled mass matrix, we are now able to
# define our problem in mass matrix form.
rhs = ODEFunction(navierstokes!, mass_matrix=M; jac_prototype=jac_sparsity)
problem = ODEProblem(rhs, u₀, (0.0,T), p);

# Now we can put everything together by specifying how to solve the problem.
# We want to use the adaptive implicit Euler method with our custom linear
# solver, which helps in the enforcement of the Dirichlet BCs. Further we
# enable the progress bar with the `progess` and `progress_steps` arguments.
# Finally we have to communicate the time step length and initialization
# algorithm. Since we start with a valid initial state we do not use one of
# DifferentialEquations.jl initialization algorithms.
# NOTE: At the time of writing this [no Hessenberg index 2 initialization is implemented](https://github.com/SciML/OrdinaryDiffEq.jl/issues/1019).
#
# To visualize the result we export the grid and our fields
# to VTK-files, which can be viewed in [ParaView](https://www.paraview.org/)
# by utilizing the corresponding pvd file.
timestepper = ImplicitEuler()
integrator = init(
    problem, timestepper, initializealg=NoInit(), dt=Δt₀,
    adaptive=true, abstol=1e-3, reltol=1e-3,
    progress=true, progress_steps=1,
    saveat=Δt_save);

pvd = paraview_collection("vortex-street.pvd");
integrator = TimeChoiceIterator(integrator, 0.0:Δt_save:T)
for (u_uc,t) in integrator
    # We ignored the Dirichlet constraints in the solution vector up to now,
    # so we have to bring them back now.
    #+
    update!(ch, t)
    u = copy(u_uc)
    apply!(u, ch)
    #compress=false flag because otherwise each vtk file will be stored in memory
    vtk_grid("vortex-street-$t.vtu", dh; compress=false) do vtk
        vtk_point_data(vtk,dh,u)
        vtk_save(vtk)
        pvd[t] = vtk
    end
end
vtk_save(pvd);

# Test the result for full proper development of the flow                   #src
using Test                                                                  #hide
function compute_divergence(dh, u, cellvalues_v)                            #hide
    divv = 0.0                                                              #hide
    @inbounds for (i,cell) in enumerate(CellIterator(dh))                   #hide
        Ferrite.reinit!(cellvalues_v, cell)                                 #hide
        for q_point in 1:getnquadpoints(cellvalues_v)                       #hide
            dΩ = getdetJdV(cellvalues_v, q_point)                           #hide
                                                                            #hide
            all_celldofs = celldofs(cell)                                   #hide
            v_celldofs = all_celldofs[dof_range(dh, :v)]                    #hide
            v_cell = u[v_celldofs]                                          #hide
                                                                            #hide
            divv += function_divergence(cellvalues_v, q_point, v_cell) * dΩ #hide
        end                                                                 #hide
    end                                                                     #hide
    return divv                                                             #hide
end                                                                         #hide
@testset "INS OrdinaryDiffEq" begin                                         #hide
    u = copy(integrator.integrator.u)                                       #hide
    apply!(u, ch)                                                           #hide
    Δdivv = abs(compute_divergence(dh, u, cellvalues_v))                    #hide
    @test isapprox(Δdivv, 0.0, atol=1e-12)                                  #hide
                                                                            #hide
    Δv = 0.0                                                                #hide
    for cell in CellIterator(dh)                                            #hide
        Ferrite.reinit!(cellvalues_v, cell)                                 #hide
        all_celldofs = celldofs(cell)                                       #hide
        v_celldofs = all_celldofs[dof_range(dh, :v)]                        #hide
        v_cell = u[v_celldofs]                                              #hide
        coords = getcoordinates(cell)                                       #hide
        for q_point in 1:getnquadpoints(cellvalues_v)                       #hide
            dΩ = getdetJdV(cellvalues_v, q_point)                           #hide
            coords_qp = spatial_coordinate(cellvalues_v, q_point, coords)   #hide
            v = function_value(cellvalues_v, q_point, v_cell)               #hide
            Δv += norm(v - parabolic_inflow_profile(coords_qp, T))^2*dΩ     #hide
        end                                                                 #hide
    end                                                                     #hide
    @test isapprox(sqrt(Δv), 0.0, atol=1e-3)                                #hide
end;                                                                        #hide

#md # ## [Plain Program](@id ns_vs_diffeq-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [ns_vs_diffeq.jl](ns_vs_diffeq.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
using CompScienceMeshes
using WiltonInts84

const p1 = point(0,0,0)
const p2 = point(1,0,0)
const p3 = point(0,1,0)

const c = point(0.2, 0.2, 0.2)
const r, R = 0.5, 1.3

function payload()
    I,K,G = wiltonints(p1,p2,p3,c,r,R,Val{1})
end

function driver()
    payload()
    Profile.clear_malloc_data()
    payload()
end

driver()

ext: jl


text:
using Test
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges

struct ObjectiveDummyBridge <: MOIB.Objective.AbstractBridge
    id::Int
end

function test_empty(map)
    @test isempty(map)
    @test length(map) == 0
    @test isempty(values(map))
    @test iterate(map) === nothing
    @test MOIB.Objective.function_type(map) === nothing
end

map = MOIB.Objective.Map()
test_empty(map)
@test sprint(MOIB.print_num_bridges, map) == "\nwith 0 objective bridges"

x = MOI.VariableIndex(1)
fx = MOI.SingleVariable(x)
MOIB.Objective.add_key_for_bridge(map, ObjectiveDummyBridge(1), fx)
@test MOIB.Objective.root_bridge(map) == ObjectiveDummyBridge(1)
@test sprint(MOIB.print_num_bridges, map) == "\nwith 1 objective bridge"
func = 1.0fx
MOIB.Objective.add_key_for_bridge(map, ObjectiveDummyBridge(2), func)
@test MOIB.Objective.root_bridge(map) == ObjectiveDummyBridge(2)
@test sprint(MOIB.print_num_bridges, map) == "\nwith 2 objective bridges"

empty!(map)
test_empty(map)
@test sprint(MOIB.print_num_bridges, map) == "\nwith 0 objective bridges"

@testset "EmptyMap" begin
    map = MOIB.Objective.EmptyMap()
    test_empty(map)
    empty!(map)
    test_empty(map)
    @test sprint(MOIB.print_num_bridges, map) == ""
end


text:
module ExtendedVariantModeExtendedVariantOnScreenDisplayTrigger

using ..Ahorn, Maple

@mapdef Trigger "ExtendedVariantMode/ExtendedVariantOnScreenDisplayTrigger" ExtendedVariantOnScreenDisplayTrigger(x::Integer, y::Integer, width::Integer=16, height::Integer=16,
	enable::Bool=true)

const placements = Ahorn.PlacementDict(
	"ExtendedVariant On-Screen Display Trigger\n(Extended Variant Mode)" => Ahorn.EntityPlacement(
		ExtendedVariantOnScreenDisplayTrigger,
		"rectangle"
	)
)

end


text:
using Random, Plots
using Zygote, ForwardDiff
using OrdinaryDiffEq, DiffEqSensitivity
using LinearAlgebra
using Statistics
using ProgressBars, Printf
using Flux
using Flux.Optimise: update!
using Flux.Losses: mae, mse
using BSON: @save, @load
using DelimitedFiles
using YAML

ENV["GKSwstype"] = "100"

cd(dirname(@__DIR__))
conf = YAML.load_file("./config.yaml")

expr_name = conf["expr_name"]
fig_path = string("./results/", expr_name, "/figs")
ckpt_path = string("./results/", expr_name, "/checkpoint")
config_path = "./results/$expr_name/config.yaml"

is_restart = Bool(conf["is_restart"])
ns = Int64(conf["ns"])
nr = Int64(conf["nr"])
lb = Float64(conf["lb"])
n_epoch = Int64(conf["n_epoch"])
n_plot = Int64(conf["n_plot"])
grad_max = Float64(conf["grad_max"])
maxiters = Int64(conf["maxiters"])

lr_max = Float64(conf["lr_max"])
lr_min = Float64(conf["lr_min"])
lr_decay = Float64(conf["lr_decay"])
lr_decay_step = Int64(conf["lr_decay_step"])
w_decay = Float64(conf["w_decay"])

llb = lb;
global p_cutoff = -1.0

const l_exp = 1:14
n_exp = length(l_exp)

l_train = []
l_val = []
for i = 1:n_exp
    j = l_exp[i]
    if !(j in [2, 6, 9, 12])
        push!(l_train, i)
    else
        push!(l_val, i)
    end
end

opt = Flux.Optimiser(
    ExpDecay(lr_max, lr_decay, length(l_train) * lr_decay_step, lr_min),
    ADAMW(0.005, (0.9, 0.999), w_decay),
);

if !is_restart
    if ispath(fig_path)
        rm(fig_path, recursive = true)
    end
    if ispath(ckpt_path)
        rm(ckpt_path, recursive = true)
    end
end

if ispath("./results") == false
    mkdir("./results")
end

if ispath("./results/$expr_name") == false
    mkdir("./results/$expr_name")
end

if ispath(fig_path) == false
    mkdir(fig_path)
    mkdir(string(fig_path, "/conditions"))
end

if ispath(ckpt_path) == false
    mkdir(ckpt_path)
end

cp("./config.yaml", config_path; force=true)


text:
export  lex_number, lex_string, lex_identifier, lex_hexbinary,
        lex_basic

"Lex a number. Could be hexadecimal, scientific or start with +-"
function lex_number(l::Lexer)
    scan_number(l)
    emit_token(l, NUMBER)
    lex_basic
end

function lex_hexbinary(l::Lexer)
    accept_char(l, "<")
    accept_char_run(l) do ch
        isxdigit(ch) || isspace(ch)
    end
    if accept_char(l, ">")
        emit_token(l, HEXBINARY, filter(isxdigit, lexeme(l)))
    else
        return error(l, "Hexadecimal data must end with >")
    end
    return lex_basic
end

"Lex a string enclosed in \" quotes"
function lex_string(l::Lexer)
    if scan_string(l)
        emit_token(l, STRING)
        lex_basic
    else
        error(l, "Not a valid quoted string")
    end
end

"lex an identifier such as a variable or function name"
function lex_identifier(l::Lexer)
    scan_identifier(l)
    emit_token(l, IDENT)
    return lex_basic
end

function lex_basic(l::Lexer)
	while true
        ignore_whitespace(l)
    	ch = peek_char(l)

    	if ch == EOFChar
    		emit_token(l, EOF)
            return lex_end
    	elseif ch in "{}(),=;"
    		next_char(l)
    		emit_token(l, TokenType(ch))
    	elseif isdigit(ch) || ch in "-+"
    		return lex_number
    	elseif ch == '"'
    		return lex_string
    	elseif isletter(ch)
    		return lex_identifier
        elseif ch == '<'
            return lex_hexbinary
    	end
    end
end


text:
#=
This page defines a number of composite types:

1. Data
    - Contains all the training, validation and test data.
    - The training data is fed in batches into the network during learning.

2. Layer
    - An abstract type with a number of concrete layer types:

    Input Layers:
    i.    InputLayer

    Hidden Layers:
    ii.   FullyConnectedLayer
    iii.  WordEmbeddingLayer
    iv.   SparseEncoderLayer
    v.    ConvolutionLayer
    vi.   MeanPoolLayer
    v.    MaxPoolLayer

    Output Layers:
    vi.   LinearOutputLayer
    vii.  MultiLinearOutputLayer
    viii. LogisticOutputLayer
    ix.   SoftmaxOutputLayer

    TODO - MultiLogistic, MultiSoftmax

3. NeuralNet
    - The complete network architecture.
    - Contains the data, network layers, and regularisation.
=#

########################################
# Data Container
########################################

abstract Data

type UData{R1<:Real} <: Data
    X_train::Array{R1}
end

function Data(X_train)
    UData(X_train)
end

type TData{R1<:Real, R2<:Real} <: Data
    X_train::Array{R1}
    y_train::Array{R2}
end

function Data(X_train, y_train)
    TData(X_train, y_train)
end

type TVData{R1<:Real, R2<:Real} <: Data
    X_train::Array{R1}
    y_train::Array{R2}
    X_val::Array{R1}
    y_val::Array{R2}
end

function Data(X_train, y_train, X_val, y_val)
    TVData(X_train, y_train, X_val, y_val)
end

type TVTData{R1<:Real, R2<:Real} <: Data
    X_train::Array{R1}
    y_train::Array{R2}
    X_val::Array{R1}
    y_val::Array{R2}
    X_test::Array{R1}
    y_test::Array{R2}
end

function Data(X_train, y_train, X_val, y_val, X_test, y_test)
    TVTData(X_train, y_train, X_val, y_val, X_test, y_test)
end

########################################
# Layer
########################################

# Define an abstract type as a parent of the layers defined below
abstract NetLayer

# Define some type unions for the input values to the layer constructors
UTupInt = Union{Tuple{Vararg{Integer}}, Integer}
USymDistn = Union{Symbol, Distribution}

########################################
# Input Layer
########################################

# Define InputLayer composite type
type InputLayer{R<:Real} <: NetLayer
    A::Array{R}       # input data
    Ar::Matrix{R}     # shared memory matrix reshape for chaining with next layer
end

# InputLayer outer constructor
function InputLayer{D<:Data, T<:Integer}(databox::D, batch_size::T)

    # Input array - convert from SubArray to Array
    X = convert(Array, viewbatch(databox.X_train, 1, batch_size))

    # Dims for reshaped 2D matrix
    rows = prod(size(X)[1:end-1])
    cols = size(X)[end]

    InputLayer(
    X, reshape(X, rows, cols)
    )
end

########################################
# Fully Connected Layer
########################################

# Define FullyConnectedLayer composite type
type FullyConnectedLayer{F<:AbstractFloat} <: NetLayer
    b::Vector{F}       # bias vector
    W::Matrix{F}       # weight matrix
    Z::Matrix{F}       # weighted sum of inputs and params i.e. Z⁽³⁾ = A⁽²⁾ W⁽³⁾ + b⁽³⁾
    A::Matrix{F}       # post activation values i.e. A⁽³⁾ = f(Z⁽³⁾)
    Δ::Matrix{F}       # delta matrix
    ∇b::Vector{F}      # bias gradient
    ∇W::Matrix{F}      # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
    f!::Function       # activation function
    ∇f!::Function      # activation function derivative and hademaard product
end

# FullyConnectedLayer outer constructor
function FullyConnectedLayer(F::DataType, input_dims::UTupInt, fan_out::Integer;
                             activation::Symbol = :logistic,
                             init::USymDistn = :default)

    # Dims
    fan_in = prod(input_dims[1:(end-1)])
    batch_size = input_dims[end]

    # Initialise weights
    distn = choose_distn(init, activation, fan_in, fan_out)
    W = F.(rand(distn, fan_out, fan_in))
    b = zeros(F, fan_out)

    # Build layer
    FullyConnectedLayer(
    b, W,
    zeros(F, fan_out, batch_size),
    zeros(F, fan_out, batch_size),
    zeros(F, fan_out, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    activation_dict[activation], activation_deriv_dict[activation]
    )
end

# Default types are 64 bit
function FullyConnectedLayer(input_dims::UTupInt, fan_out::Integer;
                             activation::Symbol = :logistic,
                             init::USymDistn = :default)
    FullyConnectedLayer(Float64, input_dims, fan_out,
                        activation = activation, init = init)
end

########################################
# Word Embedding Layer
########################################

# Define WordEmbeddingLayer composite type
type WordEmbeddingLayer{F<:AbstractFloat} <: NetLayer
    W::Matrix{F}       # weight matrix
    A::Matrix{F}       # feature values of each input n-gram
    Δ::Matrix{F}       # delta matrix
    ∇W::Matrix{F}      # weight gradient
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
end

# WordEmbeddingLayer outer constructor
function WordEmbeddingLayer(F::DataType, input_dims::UTupInt,
                            vocab_size::Integer, num_feats::Integer;
                            init::USymDistn = :default)

    # Dims
    num_words = input_dims[1]
    batch_size = input_dims[2]

    # Initialise weights
    if init == :default
        distn = Normal(0, 0.01)
    elseif init <: Symbol
        println("Error: Named initialisations are not used for Word Embedding Layer.
               Enter instead a distribution or leave out and the default
               distribution - Normal(0, 0.01) - will be used.")
        return
    elseif init <: Distribution
        distn = init
    end
    W = F.(rand(distn, num_feats, vocab_size))

    # Build layer
    WordEmbeddingLayer(
    W,
    zeros(F, num_words * num_feats, batch_size),
    zeros(F, num_words * num_feats, batch_size),
    zeros(W), zeros(W)
    )
end

# Default types are 64 bit
function WordEmbeddingLayer(input_dims::UTupInt, vocab_size::Integer,
                            num_feats::Integer; init::USymDistn = :default)
    WordEmbeddingLayer(Float64, input_dims, vocab_size, num_feats,
                       init = init)
end

########################################
# Sparse Encoder Layer
########################################

# Define SparseEncoderLayer composite type
type SparseEncoderLayer{F<:AbstractFloat} <: NetLayer
    b::Vector{F}       # bias vector
    W::Matrix{F}       # weight matrix
    Z::Matrix{F}       # weighted sum of inputs and params i.e. Z⁽³⁾ = A⁽²⁾ W⁽³⁾ + b⁽³⁾
    A::Matrix{F}       # post activation values i.e. A⁽³⁾ = f(Z⁽³⁾)
    Δ::Matrix{F}       # delta matrix
    ∇b::Vector{F}      # bias gradient
    ∇W::Matrix{F}      # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
    p̂::Vector{F}       # average activation - "prob" of each neuron firing
    spgrad::Vector{F}  # gradient of the cost of sparsity
    p::F               # sparsity parameter
    β::F               # controls the weight of the sparsity penalty term
    f!::Function       # activation function
    ∇f!::Function      # activation function derivative and hademaard product
    δ!::Function       # delta function
    m::F               # batch_size
end

# SparseEncoderLayer outer constructor
function SparseEncoderLayer(F::DataType, input_dims::UTupInt, fan_out::Integer,
                            p::Real, β::Real;
                            activation::Symbol = :logistic,
                            init::USymDistn = :default)

    # Dims
    fan_in = prod(input_dims[1:(end-1)])
    batch_size = input_dims[end]

    # Initialise weights
    distn = choose_distn(init, activation, fan_in, fan_out)
    W = F.(rand(distn, fan_out, fan_in))
    b = zeros(F, fan_out)

    # Build layer
    SparseEncoderLayer(
    b, W,
    zeros(F, fan_out, batch_size),
    zeros(F, fan_out, batch_size),
    zeros(F, fan_out, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    zeros(b), zeros(b), p, β,
    activation_dict[activation], activation_deriv_dict[activation],
    delta_sparse!,
    F(batch_size)
    )
end

# Default types are 64 bit
function SparseEncoderLayer(input_dims::UTupInt, fan_out::Integer,
                            p::Real, β::Real;
                            activation::Symbol = :logistic,
                            init::USymDistn = :default)
    SparseEncoderLayer(Float64, input_dims, fan_out, p, β,
                       activation = activation, init = init)
end

########################################
# Convolution Layer
########################################

# Define ConvolutionLayer composite type
type ConvolutionLayer{F<:AbstractFloat} <: NetLayer
    b::Vector{F}       # bias vector
    W::Array{F}        # weight array
    Z::Array{F}        # convolved sum of inputs and params plus bias i.e. Z⁽³⁾ = A⁽²⁾ ∗ W⁽³⁾ + b⁽³⁾
    Zr::Matrix{F}      # shared memory matrix reshape for chaining with next layer
    A::Array{F}        # post activation values i.e. A⁽³⁾ = f(Z⁽³⁾)
    Ar::Matrix{F}      # shared memory matrix reshape for chaining with next layer
    Δ::Array{F}        # delta array
    Δr::Matrix{F}      # shared memory matrix reshape for chaining with next layer
    ∇b::Vector{F}      # bias gradient
    ∇W::Array{F}       # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Array{F}    # weight velocity (for momentum learning)
    f!::Function       # activation function
    ∇f!::Function      # activation function derivative and hademaard product
end

# ConvolutionLayer outer constructor
function ConvolutionLayer(F::DataType, input_dims::Tuple, patch_dims::Tuple;
                          activation::Symbol = :logistic,
                          init::USymDistn = :default)

    # Dims
    fan_in = prod(patch_dims[1:(end - 1)])
    fan_out = patch_dims[end]
    conv_dims = Array(Int, 4)   # dim x dim x num_patches x num_obs
    conv_dims[1] = input_dims[1] - patch_dims[1] + 1
    conv_dims[2] = input_dims[2] - patch_dims[2] + 1
    conv_dims[3] = fan_out
    conv_dims[4] = input_dims[end]

    # Initialise weights
    distn = choose_distn(init, activation, fan_in, fan_out)
    W = F.(rand(distn, patch_dims...))
    b = zeros(F, fan_out)

    # Z, A & Δ - Initialise for shared memory reshaping
    Z = zeros(F, conv_dims...)
    A = zeros(F, conv_dims...)
    Δ = zeros(F, conv_dims...)
    rows = prod(conv_dims[1:end-1])
    cols = conv_dims[end]

    # Build layer
    ConvolutionLayer(
    b, W,
    Z, reshape(Z, rows, cols),
    A, reshape(A, rows, cols),
    Δ, reshape(Δ, rows, cols),
    zeros(b), zeros(W), zeros(b), zeros(W),
    activation_dict[activation], activation_deriv_dict[activation]
    )
end

# Default types are 64 bit
function ConvolutionLayer(input_dims::Tuple, patch_dims::Tuple;
                          activation::Symbol = :logistic,
                          init::USymDistn = :default)
    ConvolutionLayer(Float64, input_dims, patch_dims,
                     activation = activation, init = init)
end

########################################
# Mean Pooling Layer
########################################

# Define MeanPoolLayer composite type
type MeanPoolLayer{F<:AbstractFloat} <: NetLayer
    A::Array{F}        # post activation values i.e. A⁽³⁾ = meanpool(A⁽²⁾)
    Ar::Matrix{F}      # shared memory matrix reshape for chaining with next layer
    Δ::Array{F}        # delta array
    Δr::Matrix{F}      # shared memory matrix reshape for chaining with next layer
end

# MeanPoolLayer outer constructor
function MeanPoolLayer(F::DataType, input_dims::UTupInt, stride::Integer)

    # Dims
    @assert rem(input_dims[1], stride) == 0
    pool_dim = div(input_dims[1], stride)
    pool_dims = (pool_dim, pool_dim, input_dims[end - 1], input_dims[end])

    # A & Δ - Initialise for shared memory reshaping
    A = zeros(F, pool_dims)
    Δ = zeros(F, pool_dims)
    rows = prod(pool_dims[1:(end - 1)])
    cols = pool_dims[end]

    # Build layer
    MeanPoolLayer(
    A, reshape(A, rows, cols),
    Δ, reshape(Δ, rows, cols)
    )
end

# Default types are 64 bit
function MeanPoolLayer(input_dims::UTupInt, stride::Integer)
    MeanPoolLayer(Float64, input_dims, stride)
end

########################################
# Max Pooling Layer
########################################

# Define MaxPoolLayer composite type
type MaxPoolLayer{F<:AbstractFloat, T<:Integer} <: NetLayer
    A::Array{F}        # post activation values i.e. A⁽³⁾ = maxpool(A⁽²⁾)
    Ar::Matrix{F}      # shared memory matrix reshape for chaining with next layer
    Δ::Array{F}        # delta array
    Δr::Matrix{F}      # shared memory matrix reshape for chaining with next layer
    maxinds::Array{T}  # linear index of the max value from each patch
end

# MaxPoolLayer outer constructor
function MaxPoolLayer(F::DataType, input_dims::UTupInt, stride::Integer)

    # Dims
    @assert rem(input_dims[1], stride) == 0
    pool_dim = div(input_dims[1], stride)
    pool_dims = (pool_dim, pool_dim, input_dims[end - 1], input_dims[end])

    # A & Δ - Initialise for shared memory reshaping
    A = zeros(F, pool_dims)
    Δ = zeros(F, pool_dims)
    rows = prod(pool_dims[1:end-1])
    cols = pool_dims[end]

    # Build layer
    MaxPoolLayer(
    A, reshape(A, rows, cols),
    Δ, reshape(Δ, rows, cols),
    zeros(Int, input_dims...)
    )
end

# Default types are 64 bit
function MaxPoolLayer(input_dims::UTupInt, stride::Integer)
    MaxPoolLayer(Float64, input_dims, stride)
end

########################################
# Linear Output Layer
########################################

# Define LinearOutputLayer composite type
type LinearOutputLayer{F<:AbstractFloat} <: NetLayer
    b::Vector{F}       # bias vector
    W::Matrix{F}       # weight matrix
    A::Matrix{F}       # weighted sum of inputs and weights plus bias - Z⁽ᴸ⁾ = A⁽ᴸ⁻¹⁾ W⁽ᴸ⁾ + b⁽ᴸ⁾
    Δ::Matrix{F}       # delta matrix
    ∇b::Vector{F}      # bias gradient
    ∇W::Matrix{F}      # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
    δ!::Function       # delta function
    y::Vector{F}       # target vector / labels
    cost::Function     # cost function
    m::F               # batch_size
end

# LinearOutputLayer outer constructor
function LinearOutputLayer(F::DataType, databox::Data, input_dims::UTupInt;
                           init::Distribution = Bernoulli(0))

    # Dims
    dim_in = prod(input_dims[1:(end-1)])
    batch_size = input_dims[end]

    # Target vector - convert from SubArray to Array
    y = convert(Array, viewbatch(databox.y_train, 1, batch_size))

    # Check target is a vector of reals
    @assert typeof(y) <: Vector
    @assert eltype(y) <: AbstractFloat

    # Initialise weights
    W = F.(rand(init, 1, dim_in))
    b = zeros(F, 1)

    # Build layer
    LinearOutputLayer(
    b, W,
    zeros(F, 1, batch_size), zeros(F, 1, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    delta_linear!, y, quadratic_cost,
    F(batch_size)
    )
end

# Default types are 64 bit
function LinearOutputLayer(databox::Data, input_dims::UTupInt;
                           init::Distribution = Bernoulli(0))
    LinearOutputLayer(Float64, databox, input_dims, init = init)
end

########################################
# Multi Linear Output Layer
########################################

# Define MultiLinearOutputLayer composite type
type MultiLinearOutputLayer{F<:AbstractFloat} <: NetLayer
    b::Vector{F}       # bias vector
    W::Matrix{F}       # weight matrix
    A::Matrix{F}       # weighted sum of inputs and weights plus bias - Z⁽ᴸ⁾ = A⁽ᴸ⁻¹⁾ W⁽ᴸ⁾ + b⁽ᴸ⁾
    Δ::Matrix{F}       # delta matrix
    ∇b::Vector{F}      # bias gradient
    ∇W::Matrix{F}      # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
    δ!::Function       # delta function
    y::Matrix{F}       # target matrix - ℝ num_outputs x num_obs
    cost::Function     # cost function
    m::F               # batch_size
end

# MultiLinearOutputLayer outer constructor
function MultiLinearOutputLayer(F::DataType, databox::Data, input_dims::UTupInt;
                                init::Distribution = Bernoulli(0))

    # Dims
    dim_in = prod(input_dims[1:(end-1)])
    batch_size = input_dims[end]

    # Target vector - convert from SubArray to Array
    y = convert(Array, viewbatch(databox.y_train, 1, batch_size))

    # Check target is a matrix of reals
    @assert typeof(y) <: Matrix
    @assert eltype(y) <: AbstractFloat

    # Number of Outputs
    num_outputs = size(y, 1)

    # Initialise weights
    W = F.(rand(init, num_outputs, dim_in))
    b = zeros(F, num_outputs)

    # Build layer
    MultiLinearOutputLayer(
    b, W,
    zeros(F, num_outputs, batch_size), zeros(F, num_outputs, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    delta_linear!, y, quadratic_cost,
    F(batch_size)
    )
end

# Default types are 64 bit
function MultiLinearOutputLayer(databox::Data, input_dims::UTupInt;
                                init::Distribution = Bernoulli(0))
    MultiLinearOutputLayer(Float64, databox, input_dims, init = init)
end

########################################
# Logistic Output Layer
########################################

# Define LogisticOutputLayer composite type
type LogisticOutputLayer{F<:AbstractFloat, T<:Integer} <: NetLayer
    b::Vector{F}       # bias vector
    W::Matrix{F}       # weight matrix
    Z::Matrix{F}       # weighted sum of inputs and weights plus bias i.e. Z⁽ᴸ⁾ = A⁽ᴸ⁻¹⁾ W⁽ᴸ⁾ + b⁽ᴸ⁾
    A::Matrix{F}       # post activation values i.e. A⁽ᴸ⁾ = f(Z⁽ᴸ⁾)
    Δ::Matrix{F}       # delta matrix
    ∇b::Vector{F}      # bias gradient
    ∇W::Matrix{F}      # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
    f!::Function       # hypothesis function
    δ!::Function       # delta function
    y::Vector{T}       # target vector / labels
    cost::Function     # cost function
    m::F               # batch_size
end

# OutputLayer outer constructor
function LogisticOutputLayer(F::DataType, databox::Data, input_dims::UTupInt;
                             init::Distribution = Bernoulli(0))

    # Dims
    dim_in = prod(input_dims[1:(end-1)])
    batch_size = input_dims[end]

    # Target vector - convert from SubArray to Array
    y = convert(Array, viewbatch(databox.y_train, 1, batch_size))

    # Check target is a vector of 0 and 1 integers
    @assert typeof(y) <: Vector
    @assert eltype(y) <: Integer
    @assert sort(unique(y)) == [0, 1]

    # Initialise weights
    W = F.(rand(init, 1, dim_in))
    b = zeros(F, 1)

    # Build layer
    LogisticOutputLayer(
    b, W,
    zeros(F, 1, batch_size), zeros(F, 1, batch_size), zeros(F, 1, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    logistic!, delta_logistic!, y, cross_entropy_cost,
    F(batch_size)
    )
end

# Default types are 64 bit
function LogisticOutputLayer(databox::Data, input_dims::UTupInt;
                             init::Distribution = Bernoulli(0))
    LogisticOutputLayer(Float64, databox, input_dims, init = init)
end

########################################
# Softmax Output Layer
########################################

# Define SoftmaxOutputLayer composite type
type SoftmaxOutputLayer{F<:AbstractFloat, T<:Integer} <: NetLayer
    b::Vector{F}       # bias vector
    W::Matrix{F}       # weight matrix
    Z::Matrix{F}       # weighted sum of inputs and weights plus bias i.e. Z⁽ᴸ⁾ = A⁽ᴸ⁻¹⁾ W⁽ᴸ⁾ + b⁽ᴸ⁾
    u::Matrix{F}       # Max row vector used in softmax function
    A::Matrix{F}       # post activation values i.e. A⁽ᴸ⁾ = f(Z⁽ᴸ⁾)
    Δ::Matrix{F}       # delta matrix
    ∇b::Vector{F}      # bias gradient
    ∇W::Matrix{F}      # weight gradient
    b_vel::Vector{F}   # bias velocity (for momentum learning)
    W_vel::Matrix{F}   # weight velocity (for momentum learning)
    f!::Function       # hypothesis function
    δ!::Function       # delta function
    y::Vector{T}       # target vector / labels
    cost::Function     # cost function
    m::F               # batch_size
end

# SoftmaxOutputLayer outer constructor
function SoftmaxOutputLayer(F::DataType, databox::Data, input_dims::UTupInt,
                            num_classes::Integer;
                            init::Distribution = Bernoulli(0))

    # Dims
    dim_in = prod(input_dims[1:(end-1)])
    batch_size = input_dims[end]

    # Target vector - convert from SubArray to Array
    y = convert(Array, viewbatch(databox.y_train, 1, batch_size))

    # Check target is a vector of integers
    @assert typeof(y) <: Vector
    @assert eltype(y) <: Integer

    # Initialise weights
    W = F.(rand(init, num_classes, dim_in))
    b = zeros(F, num_classes)

    # Build layer
    SoftmaxOutputLayer(
    b, W,
    zeros(F, num_classes, batch_size), zeros(F, 1, batch_size),
    zeros(F, num_classes, batch_size), zeros(F, num_classes, batch_size),
    zeros(b), zeros(W), zeros(b), zeros(W),
    softmax!, delta_softmax!, y, log_loss_cost,
    F(batch_size)
    )
end

# Default types are 64 bit
function SoftmaxOutputLayer(databox::Data, input_dims::UTupInt,
                            num_classes::Integer;
                            init::Distribution = Bernoulli(0))
    SoftmaxOutputLayer(Float64, databox, input_dims, num_classes, init = init)
end

########################################
# Layer Type Unions
########################################

# Output layers
LinOutputLayer = Union{LinearOutputLayer, MultiLinearOutputLayer}
OutputLayer = Union{LinOutputLayer, LogisticOutputLayer, SoftmaxOutputLayer}

# Pooling layers
PoolLayer = Union{MeanPoolLayer, MaxPoolLayer}

# High dimensional layers
HDL = Union{InputLayer, ConvolutionLayer, PoolLayer}

# Low dimensional layers
LDL = Union{FullyConnectedLayer, WordEmbeddingLayer, SparseEncoderLayer,
            OutputLayer}

########################################
# Whole Neural Network
########################################

# Define Network composite type that contains the layers
type NeuralNet{D<:Data, L<:NetLayer, F<:AbstractFloat, T<:Integer}
    dtype::DataType
    data::D
    layers::Vector{L}
    reg_cost::Function
    reg_grad!::Function
    λ::F
    batch_size::T
    last_train_loss::Vector{F}
    full_train_loss::Vector{F}
    val_loss::Vector{F}
end

# NeuralNet outer constructor
function NeuralNet{D<:Data, L<:NetLayer}(data::D, layers::Vector{L},
                                         λ=0.0; regularisation = :L2)

    # Check that the layer dimensions allow chaining
    checkdims(layers)

    # Check floating point data type is consistent
    F = eltype(layers[end].A)
    for ℓ in layers[2:end]
        @assert F == eltype(ℓ.A)
    end

    # Create Neural Network
    NeuralNet(
    F, data, layers,
    reg_cost_dict[regularisation], reg_grad_dict[regularisation], F(λ),
    size(layers[1])[end],
    Array(F, 0), Array(F, 0), Array(F, 0)
    )
end

########################################
# Check dimensions before building network
########################################

# Low -> Low
function checkdims(ℓ::LDL, prev_ℓ::LDL)
    @assert size(ℓ.W, 2) == size(prev_ℓ.A, 1)  # Number of features match weight matrix
    @assert size(ℓ.A, 2) == size(prev_ℓ.A, 2)  # Number of observations / batch size
end

# High -> Low
function checkdims(ℓ::LDL, prev_ℓ::HDL)
    @assert size(ℓ.W, 2) == size(prev_ℓ.Ar, 1)  # Number of features match weight matrix
    @assert size(ℓ.A, 2) == size(prev_ℓ.Ar, 2)  # Number of observations / batch size
end

# High -> High:

# High -> Word Embedding
function checkdims(ℓ::WordEmbeddingLayer, prev_ℓ::HDL)
    @assert size(ℓ.A, 1) / size(ℓ.W, 1) == size(prev_ℓ.A, 1)  # Number of words chains with weight matrix
    @assert size(ℓ.A)[end] == size(prev_ℓ.A)[end]  # Number of observations / batch size
end

# High -> Convolution
function checkdims(ℓ::ConvolutionLayer, prev_ℓ::HDL)
    @assert size(ℓ.A)[end] == size(prev_ℓ.A)[end]  # Number of observations / batch size
end

# High -> Pool
function checkdims(ℓ::PoolLayer, prev_ℓ::HDL)
    @assert size(ℓ.A)[end] == size(prev_ℓ.A)[end]  # Number of observations / batch size
    @assert size(prev_ℓ.A, 1) % size(ℓ.A, 1) == 0  # Pooling divides previous layer rows evenly
    @assert size(prev_ℓ.A, 2) % size(ℓ.A, 2) == 0  # Pooling divides previous layer cols evenly
end

# All
function checkdims(layers::Vector{NetLayer})
    for (ℓ, prev_ℓ) in zip(layers[2:end], layers[1:end-1])
        checkdims(ℓ, prev_ℓ)
    end
end


text:
function prod_dim{T, N}(x::Array{T, N})
s = 1
for i = 1:N
s = s * size(x, i)
end
return s
end

max_stars_repo_path: Module 2/Chapter04/prod_tim.jl

max_issues_repo_path: Module 2/Chapter04/prod_tim.jl

max_forks_repo_path: Module 2/Chapter04/prod_tim.jl


text:
function MQKP(capacity::Array{Int64,1},
              Q::Array{Int64,2},
              abort::Int64,
              solv::String,
              show_opt::Int64,
              cpu_cores::Int64,
              allowed_gap::Float64,
              max_nodes::Int64)

    if CHECKCAPACITY(Q,capacity) == 1
        mqkp = Model(GAMS.Optimizer)
        if show_opt == 0
            set_silent(mqkp)
        end
        set_optimizer_attribute(mqkp, GAMS.ModelType(), "MIQCP")
        set_optimizer_attribute(mqkp, "Solver",  solv)
        set_optimizer_attribute(mqkp, "OptCR",   allowed_gap)
        set_optimizer_attribute(mqkp, "ResLim",  abort)
        set_optimizer_attribute(mqkp, "Threads", cpu_cores)
        set_optimizer_attribute(mqkp, "NodLim",  max_nodes)
        GI = 1:size(Q,2)
        GK = 1:size(capacity,1)
        @variable(mqkp, X[GI,GK], Bin)
        @objective(mqkp, Max, sum(X[i,k] .* X[j,k] .* Q[i,j] for i in GI, j in 1:i-1, k in GK))
        @constraint(mqkp, capconstraint[k in GK], sum(X[i,k] for i in GI) == capacity[k])
        @constraint(mqkp, minallocation[i in GI], sum(X[i,k] for k in GK) >= 1)
        JuMP.optimize!(mqkp)
        G = abs(objective_bound(mqkp)-objective_value(mqkp))/abs(objective_value(mqkp)+0.00000000001)
        out = Array{Int64,2}(undef,size(Q,2),size(capacity,1)) .= 0
        for i = 1:size(out,1)
            for j = 1:size(out,2)
                if value.(X[i,j]) > 0
                    out[i,j] = 1
                end
            end
        end
        return out::Array{Int64,2},G::Float64
    end
end

max_stars_repo_name: beyondatlas/splitdeliveries

max_issues_repo_name: beyondatlas/splitdeliveries

max_forks_repo_name: beyondatlas/splitdeliveries


text: between(a, b) = collect(a:b)

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondHamiltonianSum <: AbstractBondHamiltonian{L,N}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - generator function
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct BondHamiltonianSum{L,N,H1<:AbstractBondHamiltonian{L,N},H2<:AbstractBondHamiltonian{L,N}} <: AbstractBondHamiltonian{L,N}

    # Hamiltonian 1
    h1 :: H1

    # Hamiltonian 2
    h2 :: H2

end

# export the type
export BondHamiltonianSum

################################################################################
#
#	INTERFACING / ACCESSING BOND HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# get bond term of Hamiltonian
function bondterm(
            h :: BondHamiltonianSum{L,N,H1,H2},
            b :: AbstractBond{L,NB}
        ) :: Matrix{Complex{Float64}} where {L,N,NB,H1,H2}

    # return the sum of the two sub-hamiltonians
    return bondterm(h.h1,b) .+ bondterm(h.h2,b)
end

# export the function
export bondterm

# Overwrite the base.+ method

# FALLBACK 1 - no agreeing (spin) dimension
function Base.:+(
            h1 :: H1,
            h2 :: H2
        ) :: BondHamiltonianSum{L1,N1,H1,H2} where {L1,L2,N1,N2,H1<:AbstractBondHamiltonian{L1,N1},H2<:AbstractBondHamiltonian{L2,N2}}

    # throw an error
    error("Hamiltonians 1 and 2 don't agree in (spin) dimension: N1="*string(N1)*", N2="*string(N2))
end

# FALLBACK 2 - no agreeing labels
function Base.:+(
            h1 :: H1,
            h2 :: H2
        ) :: BondHamiltonianSum{L1,N,H1,H2} where {L1,L2,N,H1<:AbstractBondHamiltonian{L1,N},H2<:AbstractBondHamiltonian{L2,N}}

    # throw an error
    error("Hamiltonians 1 and 2 don't work on the same label type: L1="*string(L1)*", L2="*string(L2))
end

# Fitting case of same L and N
function Base.:+(
            h1 :: H1,
            h2 :: H2
        ) :: BondHamiltonianSum{L,N,H1,H2} where {L,N,H1<:AbstractBondHamiltonian{L,N},H2<:AbstractBondHamiltonian{L,N}}

    # return a new bond hamiltonian
    return BondHamiltonianSum{L,N,H1,H2}(h1,h2)
end

# CUSTOM SHOWING

# all labels
function show(io::IO, h::H) where {L,N,H<:BondHamiltonianSum{L,N}}
    print(io, "Sum (Bond) Hamiltonian for two Hamiltonians which are:\n1) ")
    show(io, h.h1)
    print("\n2) ")
    show(io, h.h2)
end


text:
cd("/home/astupidbear/.julia/v0.4/JuLearn/src/ArrayFire/examples/mlp/classification/")
ENV["JuLearnBackend"]="ArrayFire"
using JuLearn,ArrayFire#,Plots
@Eval lr=1e-2 momentum=0. batchsize=50 λ=0. ▽max=0.2
include("dataProvider.jl")
include("predict.jl")

inputBase,outputBase=dataProvider();
net=chain([:d,:h,:h,:p,:l],[2 2;2 10;10 2;2 2;2 2],:classification);
@time err=backprop(net,AFArray(inputBase),AFArray(outputBase),nIter=1000)

using PyPlot;plot(err)


text:
using Jecco
using Jecco.AdS5_3_1
using Plots
gr()
include("../src/AdS5_3_1/GW_2.jl")
#/Users/apple/Dropbox/CollisionNewPotential/3+1 project/runs/spinodal2D_model1_en0.9_run01/#
dirname = "/Users/apple/Documents/Jecco.jl/data/spinodal2D_e_0.817_Lxy_30_adaptive"
field = :energy
dt = 0.1
alg = AB3()
h, h_t, Tk, kx, ky = solve_GW(dirname, field, dt=dt, alg=alg)

#=
tmax = (length(h[:,1,1])-1)*dt
t = 0:dt:tmax
imax = 1:length(t);
j=15;
plot(t[imax], real.(h[imax,j,1]), label="h", lw=3)
plot!(t[imax], real.(h_t[imax,j,1]), label="h_t", lw=3)
plot!(t[imax], real.(Tk[imax,j,1]), label="Tk", lw=3)
=#


text:
using DiffEqSensitivity, Test

A = reshape([1,0,2,3],2,2)
function f_morris(p)
    A*p
end

function linear_batch(X)
    A= 7
    B= 0.1
    @. A*X[1,:]+B*X[2,:]
end

function ishi(X)
    A= 7
    B= 0.1
    sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1])
end

function ishi_linear(X)
    A= 7
    B= 0.1
    [sin(X[1]) + A*sin(X[2])^2+ B*X[3]^4 *sin(X[1]),A*X[1]+B*X[2]]
end

function ishi_linear_batch(X)
    A= 7
    B= 0.1
    X1 = @. sin(X[1,:]) + A*sin(X[2,:])^2+ B*X[3,:]^4 *sin(X[1,:])
    X2 = @. A*X[1,:]+B*X[2,:]
    vcat(X1',X2')
end

lb = -ones(4)*π
ub = ones(4)*π

m = gsa(f_morris,Morris(p_steps=[10,10],total_num_trajectory=1000,num_trajectory=150),[[1,5],[1,5]])
@test m.means[:,1] ≈ A[:,1] atol=1e-12
@test m.means[:,2] ≈ A[:,2] atol=1e-12
@test m.variances ≈ reshape([0,0,0,0],2,2) atol=1e-12

m = gsa(f_morris,Morris(p_steps=[10,10],relative_scale=true,total_num_trajectory=1000,num_trajectory=150),[[1,5],[1,5]])
@test m.means[2,1] ≈ 0 atol=1e-12
@test m.variances[2,1] ≈ 0 atol=1e-12
@test m.means[1,2] < m.means[2,2]

m = gsa(linear_batch,Morris(p_steps=[10,10],num_trajectory=10000),[[1,5],[1,5]],batch=true)
@test m.means ≈ [7.0  0.1] atol = 1e-2

m = gsa(ishi, Morris(num_trajectory=100000), [[lb[i],ub[i]] for i in 1:4])
@test m.means[1,:] ≈ [2.25341,4.40246,2.5049,0.0] atol = 5e-2

m = gsa(ishi_linear,Morris(num_trajectory=100000),[[lb[i],ub[i]] for i in 1:4])
@test m.means[1,:] ≈ [2.25341,4.40246,2.5049,0.0] atol = 5e-2
@test m.means[2,:] ≈ [7.0,0.1,0.0,0.0] rtol = 1e-12

m = gsa(ishi_linear_batch,Morris(num_trajectory=100000),[[lb[i],ub[i]] for i in 1:4],batch=true)
@test m.means[1,:] ≈ [2.25341,4.40246,2.5049,0.0] atol = 5e-2
@test m.means[2,:] ≈ [7.0,0.1,0.0,0.0] rtol = 1e-12


text:
using Polymake

function sumAll2(NPS)
    # This computes the Minkowski sum of an array of polytopes.
    vtxmtx = vertexRepresentation(NPS[1])
    for i = 2:length(NPS)
        vtxmtx2 = vertexRepresentation(NPS[i])
        C = [
            vtxmtx[i, :] + vtxmtx2[j, :] for i = 1:size(vtxmtx, 1)
            for j = 1:size(vtxmtx2, 1)
        ]
        vtxsum = hcat(ones(length(C)), vcat(C'...))
        global P = polytope.Polytope(POINTS = vtxsum)
        vtxmtx = vertexRepresentation(P)
        println(
            "computed the sum of " *
            string(i) *
            " out of " *
            string(length(NPS)) *
            " polytopes",
        )
    end
    P, vtxmtx
end

function newtonPolytope(f, y)
    # computes the Newton polytope of the polynomial f.
    trms = terms(f)
    varfinds = findall(ℓ -> ℓ ∈ variables(f), y)
    n = length(y)
    latPts = zeros(Int64, length(trms), n + 1)
    for i = 1:length(trms)
        v = zeros(n)
        v[varfinds] = transpose(exponents(trms[i]))
        latPts[i, :] = [1 v']
    end
    P = @pm polytope.Polytope(POINTS = latPts)
    return P
end

function getLatticePoints(P)
    # computes the lattice points in (a translated version of) P
    pts_α = P.LATTICE_POINTS_GENERATORS
    exps_α = convert(Array{Int64}, pts_α[1][:, 2:end])
    shift = abs.(minimum(vcat(exps_α, zeros(Int64, 1, size(exps_α, 2))), dims = 1))
    exps_α = exps_α .+ shift
    return exps_α
end

function vertexRepresentation(P)
    # Computes a matrix whose rows are the vertices of a lattice polytope P
    vtxmtx = P.VERTICES
    vtxmtx = convert(Array{Int64}, vtxmtx)
    vtxmtx = vtxmtx[:, 2:end]
    return vtxmtx
end

function getCodegree(A)
    # Compute the codegree of the convex hull of A
    r = 1
    @polyvar x[1:size(A, 2)]
    f = sum([prod(x .^ (A[i, :])) for i = 1:size(A, 1)])
    P = newtonPolytope(f, x)
    r = P.LATTICE_CODEGREE
    return r
end

max_stars_repo_name: simontelen/JuliaEigenvalueSolver

max_issues_repo_name: simontelen/JuliaEigenvalueSolver

max_forks_repo_name: simontelen/JuliaEigenvalueSolver


text:
module GSOcc

using DataFrames
using Infiltrator

import ..Util

# Parameter
"""
    occupations

Liste von Ämtern, die in Betracht gezogen werden.

Beispiel

["Bischof",
 "Vikar",
 "Elekt",
 "Administrator",
 "Patriarch",
 "Metropolit"]
"""
global occupations = ["Bischof", "Vikar", "Elekt", "Administrator", "Patriarch", "Metropolit"]

function regexor(lms::Array{<:AbstractString, 1})
    if length(lms) > 1
        Regex(join(lms, "|"), "i")
    else
        r".+"
    end
end

global rgxocc = regexor(occupations)

"""
    setoccupations(occ::Vector{<:AbstractString})

Setze die Liste der Ämter, die bei der Abfrage in Betracht gezogen werden.
Wenn die Liste leer ist, werden alle Ämter in Betracht gezogen.

`setoccupations()`: Gib die aktuelle Liste aus.

Beispiel

```julia
setoccupations(["Pfarrer", "Vikar"])
```
"""
function setoccupations(occ::Vector{<:AbstractString})
    global occupations = occ
    global rgxocc = regexor(occupations)
    occupations
end

setoccupations() = occupations

"""
    equivrls

Liste von Zuordnungen für als gleichwertig betrachtete Ämter
"""
equivrls = ["Gewählter Bischof" => "Elekt",
            "Erwählter Bischof" => "Elekt",
            "Fürstbischof" => "Bischof"]

"""
    setequivalentoccupations(lequivocc)

Setze Liste von Zuordnungen für als gleichwertig betrachtete Ämter.

`setequivalentoccupations()`: Gib die aktuelle Liste aus.

Beispiel

```julia
setequivalentoccupations(["Gewählter Bischof" => "Elekt",
                          "Erwählter Bischof" => "Elekt"])
```
"""
function setequivalentoccupations(lequivocc)
    global equivrls
    equivrls = [Pair(p...) for p in lequivocc]
end

setequivalentoccupations() = equivrls

const KEYS = [["ae", "ab", "ao", "at"],
              ["ae", "ab", "ao"],
              ["ae", "ao", "at"],
              ["ab", "ao", "at"],
              ["ae", "ab", "at"],
              ["ae", "ao"],
              ["ab", "ao"],
              ["ao", "at"],
              ["ae", "ab"],
              ["ae", "at"],
              ["ab", "at"],
              ["ao"],
              ["ae"],
              ["ab"],
              ["at"],
              []]

# Feldbezeichnungen in den Abfragedaten
const KEYQDPLACE = :Amtsort
const KEYQDTYPE = :Amtsart
const KEYQDBEGIN = :Amtsbeginn
const KEYQDEND = :Amtsende
const KEYQDIDMONASTERY = :ID_Kloster

# Feldbezeichnungen in der Personendatenbank
const KEYGSTYPE = "bezeichnung"
const KEYGSBEGIN = "von"
const KEYGSEND = "bis"
const KEYGSIDMONASTERY = "klosterid"

# Funktionen
function islesskey(a, b)
    findfirst(isequal(a), KEYS) > findfirst(isequal(b), KEYS)
end

"""
    evaluate!(record, dfocc::AbstractDataFrame, tolocc, occmcols)

Bewerte Daten zu Ämtern.

`dfocc`: Ämter in Abfragedaten für diese eine Person.
`tolocc`: Maximal zulässige Abweichung in Start- und Enddatum eines Amtes.
`occmcols`: Spalten in den Abfragedaten
"""
function evaluate!(record, dfocc::Union{AbstractDataFrame, Nothing}, tolocc, occmcols)
    if isnothing(dfocc) return end
    maxkey = String[]
    for row in eachrow(dfocc)
        key = evaluatesingle(record, row, tolocc, occmcols)
        if islesskey(maxkey, key)
            maxkey = key
        end
    end
    append!(record["amuster"], maxkey)
    return nothing
end

"""
    evaluatesingle!(record, row, tolocc, occmcols)

Bewerte Daten zum Amt. `row` Zeile einer Ämtertabelle.
"""
function evaluatesingle(record, row, tolocc, occmcols)::Vector{String}
    # Für die Bischöfe vor 1198 gibt es oft nur eine Angabe für das
    # Jahrhundert. "[4. Jh.]"

    matchkey = ""
    akey = String[]

    ftype = false
    fbegin = false
    fend = false
    fplace = false
    score = 0

    typeqd = KEYQDTYPE in occmcols ? row[KEYQDTYPE] : missing
    diocqd = KEYQDPLACE in occmcols ? row[KEYQDPLACE] : missing
    idmonasteryqd = KEYQDIDMONASTERY in occmcols ? row[KEYQDIDMONASTERY] : missing

    beginqd::Union{Int, Missing} = missing
    if KEYQDBEGIN in occmcols
        sdate = row[KEYQDBEGIN]
        if Util.hasdata(sdate)
            beginqd = parsedate(row[KEYQDBEGIN])
        end
    end

    endqd::Union{Int, Missing} = missing
    if KEYQDEND in occmcols
        sdate = row[KEYQDEND]
        if Util.hasdata(sdate)
            endqd = parsedate(row[KEYQDEND])
        end
    end

    maxkey = String[]
    for occrec in record["aemter"]
        key = String[]
        # Betrachte nur passende Ämter in `occupations`.
        occgs = occrec[KEYGSTYPE]
        rgm = match(rgxocc, occgs)
        rgm == nothing && continue

        fplace = false
        ftype = false
        fbegin = false
        fend = false

        # Amtsende
        if !ismissing(endqd) && endqd != "" && evaluatedate(occrec[KEYGSEND], endqd, tolocc)
            push!(key, "ae")
        end

        # Amtsbeginn
        if !ismissing(beginqd) && beginqd != "" && evaluatedate(occrec[KEYGSBEGIN], beginqd, tolocc)
            push!(key, "ab")
        end

        # Amtsort or ID_Kloster
        if (!ismissing(diocqd) && diocqd != "" && matchplace(occrec, diocqd)
            || (!ismissing(idmonasteryqd) && string(idmonasteryqd) == occrec[KEYGSIDMONASTERY]))
            push!(key, "ao")
        end

        # Amtsbezeichnung
        ftype = !ismissing(typeqd) && typeqd != "" && matchoccupation(occgs, typeqd)
        if ftype
            push!(key, "at")
        end

        if islesskey(maxkey, key)
            maxkey = key
            record["amt"] = occrec
        end
    end

    return maxkey
end

let
    global parsedate
    rgxyear = r"[0-9]?[0-9]?[0-9]{2}"

    """

    Gib `nothing` zurück, wenn kein gültiges Datum gefunden werden kann
    """
    function parsedate(sdate::Union{AbstractString, Missing})::Union{Int, Missing}
        valdate = missing
        if ismissing(sdate) return valdate end

        rgm = match(rgxyear, sdate)
        if rgm == nothing
            @warn ("Ungültiges Datum in: " * sdate)
        else
            valdate = parse(Int, rgm.match)
        end

        return valdate
    end
end

function parsedate(sdate::Int)::Int
    sdate
end

function evaluatedate(sdategs::Union{AbstractString, Missing},
                      dateqd::Union{Int, Missing},
                      tolocc)
    if ismissing(dateqd) || isnothing(dateqd) return false end
    if ismissing(sdategs) return false end
    if !Util.hasdata(sdategs) return false end
    dategs = parsedate(sdategs)
    if ismissing(dategs) return false end

    delta = abs(dategs - dateqd)
    if delta <= tolocc
        return true
    else
        return false
    end
end

"""
    matchplace(occrec, diocqd::Union{AbstractString, Missing})

Prüfe, ob die Orte in `occrec` zu der Angabe der neuen Quelle passen.
"""
function matchplace(occrec, diocqd::Union{AbstractString, Missing})
    placegs = occrec["ort"]
    diocgs = occrec["dioezese"]
    return (Util.checkname(diocqd, placegs)
            || Util.checkname(diocqd, diocgs))
end

"""
    gettypeoccupation(socc)

Extrahiere die Art des Amtes (= eines aus `occupations`)
"""
function gettypeoccupation(socc::AbstractString)
    global rgxocc
    rgm = match(rgxocc, socc)
    if rgm == nothing
        return ""
    else
        return rgm.match
    end
end

"""
    gettypeoccupation(dfocc)

Extrahiere die Art des Amtes (= eines aus `occupations`)
"""
function gettypeoccupation(dfocc::AbstractDataFrame)
    global rgxocc
    # Beginne mit dem letzen Amt
    aocc = dfocc[!, :Amtsart]
    if length(aocc) == 0
        return ""
    end
    typeocc = ""
    for occ in Iterators.reverse(aocc)
        rgm = match(rgxocc, occ)
        if rgm != nothing
            typeocc = rgm.match
            break
        end
    end
    if typeocc == ""
        typeocc = aocc[end]
    end
    typeocc
end

"""
    matchoccupation(occqd::AbstractString, occgs::AbstractString)

Vergleich Ämter direkt, oder über eine Bezeichnung für ein Bischofsamt.
"""
function matchoccupation(occgs::AbstractString, occqd::AbstractString)

    # Direkte Übereinstimmung
    Util.checkname(occqd, occgs) && return true

    global equivrls
    for rl in equivrls::Array{Pair{String, String}}
        a = replace(occgs, rl)
        b = replace(occqd, rl)
        Util.checkname(occqd, occgs) && return true
    end
    return false

end

getlastocc(record) = getlastoccinrecord(record["aemter"])

"""
    getlastoccinrecord(occs)

Gib das Amt mit dem größten Amtsende-Datum zurück
"""
function getlastoccinrecord(occs)
    maxend = 0
    occend = 0
    lastocc::Union{eltype(occs), Nothing} = nothing
    if length(occs) > 0
        lastocc = occs[1]
    end
    for occ in occs
        soccend = occ[KEYGSEND]
        if Util.hasdata(soccend)
            occend = parsedate(soccend)
            if occend > maxend
                maxend = occend
                lastocc = occ
            end
        end
    end
    lastocc
end

end # module Occ


text:
using QuandlAccess
using QuandlAccess: enhance, positive_integer, non_negative_integer, valid_date, make_url
using Dates
using Test

@testset "QuandlAccess.jl" begin

    @testset "URL construction" begin

        @test enhance("abc", "k", nothing, positive_integer) == "abc"

        @test enhance("abc", "k", 1, positive_integer) == "abc&k=1"
        @test_throws ArgumentError enhance("abc", "k", 0, positive_integer)

        @test enhance("abc", "k", 0, non_negative_integer) == "abc&k=0"
        @test_throws ArgumentError enhance("abc", "k", -1, non_negative_integer)

        @test valid_date("2020-01-01") == true
        @test valid_date("2020-1-2") == false
        @test valid_date(Date(2020,1,1)) == true
        @test valid_date(DateTime(2020,1,1,0,0,0)) == false

        let q = Quandl("abc"), s = TimeSeries("s/t"), t = Table("x/y")
            @test endswith(make_url(q, s; start_date = "2020-01-01"), "start_date=2020-01-01")
            @test endswith(make_url(q, s; end_date = "2020-01-01"), "end_date=2020-01-01")
            @test endswith(make_url(q, s; collapse = "monthly"), "collapse=monthly")
            @test endswith(make_url(q, s; transform = "diff"), "transform=diff")
            @test endswith(make_url(q, s; limit = 2), "limit=2")
            @test endswith(make_url(q, s; column_index = 2), "column_index=2")
            @test endswith(make_url(q, s; order = "asc"), "order=asc")

            @test endswith(make_url(q, t; filters=[eq("a","b")]), "a=b")
            @test endswith(make_url(q, t; filters=[lt("a","b")]), "a.lt=b")
            @test endswith(make_url(q, t; filters=[gt("a","b")]), "a.gt=b")
            @test endswith(make_url(q, t; filters=[lte("a","b")]), "a.lte=b")
            @test endswith(make_url(q, t; filters=[gte("a","b")]), "a.gte=b")
            @test endswith(make_url(q, t; columns=["a","b"]), "columns=a,b")
        end
    end

end


text:
using SparseArrays
println("testing incidence_to_triangles")
# function test_incidence_to_triangles()
D = [
    1 -1 0 0 0;
    1 0 -1 0 0;
    0 1 -1 0 0;
    0 1 0 -1 0;
    0 1 0 0 -1;
    0 0 0 1 -1
];
triangles = [
    1 2 3;
    2 4 5
];
tri = incidence_to_triangles(sparse(D'))
@test triangles' == tri

D = [
    1 -1 0 0 0;
    1 0 -1 0 0;
    1 0 0 -1 0;
    1 0 0 0 -1;
    0 1 -1 0 0;
    0 1 0 -1 0;
    0 1 0 0 -1;
    0 0 1 -1 0;
    0 0 1 0 -1;
    0 0 0 1 -1
];
triangles = [
    1 2 3;
    1 2 4;
    1 2 5;
    1 3 4;
    1 3 5;
    1 4 5;
    2 3 4;
    2 3 5;
    2 4 5;
    3 4 5
];
#tri = incidence_to_triangles(sparse(D')) #error, more than three node connect to node 1, not allowed by incidence_to_incident_list_array
#@test triangles' == tri

D = [
    1 -1 0 0 0;
    0 0 1 -1 0;
    0 0 1 0 -1;
    1 0 -1 0 0;
    1 0 0 -1 0;
    1 0 0 0 -1;
    0 1 -1 0 0;
    0 1 0 -1 0;
    0 1 0 0 -1;
    0 0 0 1 -1
];
triangles = [
    1 2 3;
    1 2 4;
    1 2 5;
    1 3 4;
    1 3 5;
    1 4 5;
    2 3 4;
    2 3 5;
    2 4 5;
    3 4 5
];
#tri = incidence_to_triangles(sparse(D')) #error, more than three node connect to node 1, not allowed by incidence_to_incident_list_array
#@test triangles' == tri

# function test_incidence_to_dict()
# D = [
#     1 -1 0 1 0;
#     0 0 1 -1 0
# ];
# node_dict = Dict(1 => Set(2), 3 => Set(4))
# nd = incidence_to_dict(D)
# @test node_dict == nd

# D = [
#     1 -1 0 1 0;
#     1 -1 1 1 0;
#     1 0 1 1 0
# ];
# node_dict = Dict(1 => Set([2, 3]))
# nd = incidence_to_dict(D)
# @test node_dict == nd


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

using Printf
using Base.Threads
import Base.Threads: @threads, nthreads, threadid

function solve!(
    mip::MIP;
    time_limit::Float64 = Inf,
    node_limit::Int = typemax(Int),
    gap_limit::Float64 = 1e-4,
    print_interval::Int = 5,
    initial_primal_bound::Float64 = Inf,
    detailed_output::Bool = false,
    branch_rule::VariableBranchingRule = HybridBranching(),
)::NodePool
    time_initial = time()
    pool = NodePool(mip = mip)
    pool.primal_bound = initial_primal_bound
    print_progress_header(detailed_output = detailed_output)

    root_node = _create_node(mip)
    if isempty(root_node.fractional_variables)
        println("root relaxation is integer feasible")
        pool.dual_bound = root_node.obj
        pool.primal_bound = root_node.obj
        return pool
    end

    offer(pool, parent_node = nothing, child_nodes = [root_node])
    @threads for t = 1:nthreads()
        child_one, child_zero, suggestions = nothing, nothing, Node[]
        while true
            time_elapsed = time() - time_initial
            if child_one !== nothing
                suggestions = Node[child_one, child_zero]
            end
            node = take(
                pool,
                suggestions = suggestions,
                time_remaining = time_limit - time_elapsed,
                node_limit = node_limit,
                gap_limit = gap_limit,
            )
            if node == :END
                break
            elseif node == :WAIT
                sleep(0.1)
                continue
            else
                ids = generate_indices(pool, 2)
                branch_variable = find_branching_var(branch_rule, node, pool)
                child_zero = _create_node(
                    mip,
                    index = ids[1],
                    parent = node,
                    branch_variable = branch_variable,
                    branch_value = 0.0,
                )
                child_one = _create_node(
                    mip,
                    index = ids[2],
                    parent = node,
                    branch_variable = branch_variable,
                    branch_value = 1.0,
                )
                offer(
                    pool,
                    parent_node = node,
                    child_nodes = [child_one, child_zero],
                    time_elapsed = time_elapsed,
                    print_interval = print_interval,
                    detailed_output = detailed_output,
                )
            end
        end
    end
    return pool
end

function _create_node(
    mip;
    index::Int = 0,
    parent::Union{Nothing,Node} = nothing,
    branch_variable::Union{Nothing,Variable} = nothing,
    branch_value::Union{Nothing,Float64} = nothing,
)::Node
    if parent === nothing
        branch_variables = Variable[]
        branch_values = Float64[]
        depth = 1
    else
        branch_variables = [parent.branch_variables; branch_variable]
        branch_values = [parent.branch_values; branch_value]
        depth = parent.depth + 1
    end
    set_bounds!(mip, branch_variables, branch_values, branch_values)
    status, obj = solve_relaxation!(mip)
    if status == :Optimal
        vals = values(mip, mip.binary_variables)
        fractional_indices =
            [j for j in 1:length(mip.binary_variables) if 1e-6 < vals[j] < 1 - 1e-6]
        fractional_values = vals[fractional_indices]
        fractional_variables = mip.binary_variables[fractional_indices]
    else
        fractional_variables = Variable[]
        fractional_values = Float64[]
    end
    n_branch = length(branch_variables)
    set_bounds!(mip, branch_variables, zeros(n_branch), ones(n_branch))
    return Node(
        mip,
        index,
        depth,
        obj,
        status,
        branch_variables,
        branch_values,
        fractional_variables,
        fractional_values,
        parent,
    )
end


text:
export BehaviorCloningPolicy

"""
    BehaviorCloningPolicy(;kw...)

# Keyword Arguments

- `approximator`: calculate the logits of possible actions directly
- `explorer=GreedyExplorer()`

"""
Base.@kwdef struct BehaviorCloningPolicy{A} <: AbstractPolicy
    approximator::A
    explorer::Any = GreedyExplorer()
end

function (p::BehaviorCloningPolicy)(env::AbstractEnv)
    s = state(env)
    s_batch = Flux.unsqueeze(s, ndims(s) + 1)
    logits = p.approximator(s_batch) |> vec  # drop dimension
    p.explorer(logits)
end

function RLBase.update!(p::BehaviorCloningPolicy, batch::NamedTuple{(:state, :action)})
    s, a = batch.state, batch.action
    m = p.approximator
    gs = gradient(params(m)) do
        ŷ = m(s)
        y = Flux.onehotbatch(a, axes(ŷ, 1))
        logitcrossentropy(ŷ, y)
    end
    update!(m, gs)
end


text:
pages = [
        "DiffEqFlux.jl: High Level Scientific Machine Learning (SciML) Pre-Built Architectures" => "index.md",
        "Differential Equation Machine Learning Tutorials" => Any[
            "examples/neural_ode.md",
            "examples/augmented_neural_ode.md",
            "examples/collocation.md",
            "examples/hamiltonian_nn.md",
            "examples/tensor_layer.md",
            "examples/multiple_shooting.md"
        ],
        "Layer APIs" => Any[
            "Classical Basis Layers" => "layers/BasisLayers.md",
            "Tensor Product Layer" => "layers/TensorLayer.md",
            "Continuous Normalizing Flows Layer" => "layers/CNFLayer.md",
            "Spline Layer" => "layers/SplineLayer.md",
            "Neural Differential Equation Layers" => "layers/NeuralDELayers.md",
            "Hamiltonian Neural Network Layer" => "layers/HamiltonianNN.md"
        ],
        "Utility Function APIs" => Any[
            "Smoothed Collocation" => "utilities/Collocation.md",
            "Multiple Shooting Functionality" => "utilities/MultipleShooting.md",
        ],
    ]


text:
# SAM File Format
# ===============

module SAM

using BioGenerics

import BioAlignments
import BioGenerics: BioGenerics, isfilled, header
import BioGenerics.Exceptions: missingerror
import BioGenerics.Automa: State
import BioSequences
import TranscodingStreams: TranscodingStreams, TranscodingStream

using Printf: @sprintf

#TODO: update import BioCore.RecordHelper: unsafe_parse_decimal
# r"[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where {T<:Unsigned}
    x = zero(T)
    @inbounds for i in range
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return x
end

# r"[-+]?[0-9]+" must match `data[range]`.
function unsafe_parse_decimal(::Type{T}, data::Vector{UInt8}, range::UnitRange{Int}) where {T<:Signed}
    lo = first(range)
    if data[lo] == UInt8('-')
        sign = T(-1)
        lo += 1
    elseif data[lo] == UInt8('+')
        sign = T(+1)
        lo += 1
    else
        sign = T(+1)
    end
    x = zero(T)
    @inbounds for i in lo:last(range)
        x = Base.Checked.checked_mul(x, 10 % T)
        x = Base.Checked.checked_add(x, (data[i] - UInt8('0')) % T)
    end
    return sign * x
end

include("flags.jl")
include("metainfo.jl")
include("record.jl")
include("header.jl")
include("reader.jl")
include("readrecord.jl")
include("writer.jl")

end


text:
# Borrowed from PyPlot.jl

###########################################################################
# Julia 0.4 help system: define a documentation object
# that lazily looks up help from a PyObject via zero or more keys.
# This saves us time when loading PyPlot, since we don't have
# to load up all of the documentation strings right away.
immutable LazyHelp
    o::PyObject
    keys::Tuple{Vararg{Compat.String}}
    LazyHelp(o::PyObject) = new(o, ())
    LazyHelp(o::PyObject, k::AbstractString) = new(o, (k,))
    LazyHelp(o::PyObject, k1::AbstractString, k2::AbstractString) = new(o, (k1,k2))
    LazyHelp(o::PyObject, k::Tuple{Vararg{AbstractString}}) = new(o, k)
end
@compat function show(io::IO, ::MIME"text/plain", h::LazyHelp)
    o = h.o
    for k in h.keys
        o = o[k]
    end
    doc = ""
    if haskey(o, "__doc__")
        ds = o["__doc__"]
        if PyCall.pynothing[] != ds.o
            doc = convert(AbstractString, ds)
        end
    end
    if isempty(doc)
        print(io, "no Python docstring found for ", o)
    else
        print(io, doc)
    end
end
Base.show(io::IO, h::LazyHelp) = @compat show(io, "text/plain", h)
function Base.Docs.catdoc(hs::LazyHelp...)
    Base.Docs.Text() do io
        for h in hs
            @compat show(io, MIME"text/plain"(), h)
        end
    end
end


text:
const SLFN = ExtremeLearningMachine.SLFN
const ActivationFunction = ExtremeLearningMachine.ActivationFunction

function prepare_model(::Type{T},
                       n_neurons::Int,
                       activation_function::ActivationFunction) where {T <: Number}
    train(T, traindata()..., n_neurons, activation_function)
end

function prepare_model(::Type{T},
                       n_samples::Int,
                       n_neurons::Int,
                       activation_function::ActivationFunction) where {T <: Number}
    train(T, traindata(n_samples)..., n_neurons, activation_function)
end

function prepare_model(::Type{T},
                       filtering_labels::Vector{Int},
                       n_neurons::Int,
                       activation_function::ActivationFunction) where {T <: Number}
    train(T, traindata(filtering_labels)..., n_neurons, activation_function)
end

function prepare_model(::Type{T},
                       n_samples::Int,
                       filtering_labels::Vector{Int},
                       n_neurons::Int,
                       activation_function::ActivationFunction) where {T <: Number}
    train(T, traindata(n_samples, filtering_labels)..., n_neurons, activation_function)
end

function train(::Type{T1},
               samples::T2,
               targets::T3,
               n_neurons::Int,
               activation_function::ActivationFunction) where {T1 <: Number, T2, T3 <: AbstractMatrix}
    d, N = size(samples)
    q, _ = size(targets)
    elm = ELM{T1}(d, q, n_neurons, activation_function)
    add_data!(elm, samples, targets)
    solve(elm)
end

function calculate_error(model::SLFN, samples::T1, truth::T2) where {T1, T2 <: AbstractMatrix}
    X = samples
    T = truth
    Y = predict(model, X)
    d, N = size(X)

    error = 1 - (sum([argmax(T[:,n]) == argmax(Y[:,n]) for n ∈ 1:N]) / N)
end


text:
using ProximalOperators
using BenchmarkTools
using LinearAlgebra
using SparseArrays
using Random

const SUITE = BenchmarkGroup()

k = "IndPSD"
SUITE[k] = BenchmarkGroup(["IndPSD"])
for T in [Float64, Complex{Float64}]
    for n in [10, 20, 50]
        SUITE[k][T, n] = @benchmarkable prox!(Y, f, X) setup=begin
            f = IndPSD()
            W = if $T <: Real
                Symmetric
            else
                Hermitian
            end
            Random.seed!(0)
            A = randn($T, $n, $n)
            X = W((A + A')/2)
            Y = similar(X)
        end
    end
end

k = "IndBox"
SUITE[k] = BenchmarkGroup(["IndBox"])
for T in [Float32, Float64]
    SUITE[k][T] = @benchmarkable prox!(y, f, x) setup=begin
        low = -ones($T, 10000)
        upp = +ones($T, 10000)
        f = IndBox(low, upp)
        x = [-2*ones($T, 3000); zeros($T, 4000); ones($T, 3000)]
        y = similar(x)
    end
end

k = "IndNonnegative"
SUITE[k] = BenchmarkGroup(["IndNonnegative"])
for T in [Float32, Float64]
    SUITE[k][T] = @benchmarkable prox!(y, f, x) setup=begin
        f = IndNonnegative()
        x = [-2*ones($T, 3000); zeros($T, 4000); ones($T, 3000)]
        y = similar(x)
    end
end

k = "NormL2"
SUITE[k] = BenchmarkGroup(["NormL2"])
for T in [Float32, Float64]
    SUITE[k][T] = @benchmarkable prox!(y, f, x) setup=begin
        f = NormL2()
        x = ones($T, 10000)
        y = similar(x)
    end
end

k = "LeastSquares"
SUITE[k] = BenchmarkGroup(["LeastSquares"])
for (T, s, sparse, iterative) in Iterators.product(
    [Float64, ComplexF64],
    [(5, 11), (11, 5)],
    [false, true],
    [false, true],
)
    SUITE[k][(T, s, sparse, iterative)] = @benchmarkable prox!(y, f, x) setup=begin
        A = if $sparse sparse(ones($T, $s)) else ones($T, $s) end
        b = ones($T, $(s[1]))
        f = LeastSquares(A, b, iterative=$iterative)
        x = ones($T, $(s[2]))
        y = similar(x)
    end
end

k = "IndExpPrimal"
SUITE[k] = BenchmarkGroup(["IndExpPrimal"])
for T in [Float32, Float64]
    SUITE[k][T] = @benchmarkable prox!(y, f, x) setup=begin
        f = IndExpPrimal()
        x = [0.537667139546100, 1.833885014595086, -2.258846861003648]
        y = similar(x)
    end
end

k = "IndSimplex"
SUITE[k] = BenchmarkGroup(["IndSimplex"])
for T in [Float32, Float64]
    SUITE[k][T] = @benchmarkable prox!(y, f, x) setup=begin
        f = IndSimplex()
        x = collect($T, -0.5:0.001:2.0)
        y = similar(x)
    end
end

k = "IndBallL1"
SUITE[k] = BenchmarkGroup(["IndBallL1"])
for T in [Float32, Float64]
    SUITE[k][T] = @benchmarkable prox!(y, f, x) setup=begin
        f = IndBallL1()
        x = collect($T, -2.0:0.001:0.5)
        y = similar(x)
    end
end


text:
# Generate aliases for all defined time scales so we can use
# e.g. `TTEpoch` instead of `Epoch{TT}`
for (scale, acronym) in zip(TimeScales.NAMES, TimeScales.ACRONYMS)
    epoch = Symbol(acronym, "Epoch")
    name = string(epoch)
    @eval begin
        const $epoch = Epoch{$scale}
        export $epoch

        """
            $($name)(str[, format])

        Construct a $($name) from a string `str`. Optionally a `format` definition can be
        passed as a [`DateFormat`](https://docs.julialang.org/en/stable/stdlib/Dates/#Dates.DateFormat)
        object or as a string. In addition to the character codes supported by `DateFormat`
        the code `D` is supported which is parsed as "day of year" (see the example below).
        The default format is `yyyy-mm-ddTHH:MM:SS.sss`.

        ### Example ###

        ```jldoctest; setup = :(using AstroTime)
        julia> $($name)("2018-02-06T20:45:00.0")
        2018-02-06T20:45:00.000 $($acronym)

        julia> $($name)("February 6, 2018", "U d, y")
        2018-02-06T00:00:00.000 $($acronym)

        julia> $($name)("2018-37T00:00", "yyyy-DDDTHH:MM")
        2018-02-06T00:00:00.000 $($acronym)
        ```
        """
        $epoch(::AbstractString)

        """
            $($name)(jd1::T, jd2::T=zero(T); origin=:j2000) where T<:AstroPeriod

        Construct a $($name) from a Julian date (optionally split into
        `jd1` and `jd2`). `origin` determines the variant of Julian
        date that is used. Possible values are:

        - `:j2000`: J2000 Julian date, starts at `2000-01-01T12:00`
        - `:julian`: Julian date, starts at `-4712-01-01T12:00`
        - `:modified_julian`: Modified Julian date, starts at `1858-11-17T00:00`

        ### Examples ###

        ```jldoctest; setup = :(using AstroTime)
        julia> $($name)(0.0days, 0.5days)
        2000-01-02T00:00:00.000 $($acronym)

        julia> $($name)(2.451545e6days, origin=:julian)
        2000-01-01T12:00:00.000 $($acronym)
        ```
        """
        $epoch(::Number, ::Number)

        """
            $($name)(year, month, day, hour=0, minute=0, second=0.0)

        Construct a $($name) from date and time components.

        ### Example ###

        ```jldoctest; setup = :(using AstroTime)
        julia> $($name)(2018, 2, 6, 20, 45, 0.0)
        2018-02-06T20:45:00.000 $($acronym)

        julia> $($name)(2018, 2, 6)
        2018-02-06T00:00:00.000 $($acronym)
        ```
        """
        $epoch(::Int, ::Int, ::Int)
    end
end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
#####
# ListNode
#####
import LinearAlgebra
Base.@kwdef mutable struct ListNode{V}
    val::V = 0
    next::Union{Nothing,ListNode{V}} = nothing
end

ListNode(x) = ListNode(x, nothing)

function ListNode{V}(x::Vector{V}) where {V}
    head = cur = ListNode(x[1])
    for i in 2:length(x)
        cur.next = ListNode(x[i])
        cur = cur.next
    end
    return head
end

val(x::ListNode) = x.val
val!(x::ListNode, y) = x.val = y
next(x::ListNode) = x.next
next!(x::ListNode, y) = x.next = y

Base.:(==)(x::ListNode, y::ListNode) = x.val == y.val && x.next == y.next

#####
# TreeNode
#####

Base.@kwdef mutable struct TreeNode{V}
    val::V = 0
    left::Union{Nothing,TreeNode} = nothing
    right::Union{Nothing,TreeNode} = nothing
end

TreeNode(x) = TreeNode(; val=x)
TreeNode(x, left) = TreeNode(; val=x, left=left)

function Base.:(==)(t1::TreeNode, t2::TreeNode)
    return t1.val == t2.val && t1.left == t2.left && t1.right == t2.right
end

function Base.convert(::Type{TreeNode{V}}, xs::Vector) where {V}
    @assert length(xs) > 0
    root = _build_tree!(V, xs)
    return root
end

TreeNode{V}(xs::Vector) where {V} = convert(TreeNode{V}, xs)

function _build_tree!(::Type{T}, xs::Vector) where T
    res = TreeNode(xs[1])
    q1 = TreeNode{T}[res]
    q2 = TreeNode{T}[]
    idx, len = 2, length(xs)
    while idx <= len
        while !isempty(q1)
            qf = popfirst!(q1)
            idx > len && break
            if isnothing(xs[idx])
                qf.left = nothing
            else
                qf.left = TreeNode(xs[idx])
                push!(q2, qf.left)
            end
            (idx += 1) > len && break
            if isnothing(xs[idx])
                qf.right = nothing
            else
                qf.right = TreeNode(xs[idx])
                push!(q2, qf.right)
            end
            idx += 1
        end
        q1, q2 = q2, q1
    end
    return res
end

function next_perm!(itr)::Bool
    (isempty(itr) || length(itr) == 1) && return false

    i = findlast(idx -> itr[idx + 1] > itr[idx], 1:lastindex(itr) - 1)
    isnothing(i) && (reverse!(itr); return false)

    j = findlast(idx -> itr[idx] > itr[i], i+1:lastindex(itr)) + i
    itr[i], itr[j] = itr[j], itr[i]
    reverse!(itr, i + 1)
    return true
end

"""
``res = mat ^ p \\pmod m``
"""
function mat_fast_mul(mat::AbstractMatrix{T}, p::T, m::T=typemax(T)) where T
    p == 0 && return (ones(T, size(mat)))
    (m == 1 || m == -1) && return (zeros(T, size(mat)))
    base = mat .% m
    res = copy(base)
    tmp = typeof(mat)(undef, size(mat))
    t = prevpow(2, p)
    p -= t
    while true
        if p >= t
            LinearAlgebra.mul!(tmp, res, base)
            tmp, res = res, tmp
            res .%= m
            p -= t
        end
        t >>>= 1
        t <= 0 && break
        LinearAlgebra.mul!(tmp, res, res)
        tmp, res = res, tmp
        res .%= m
    end
    res
end