text:
# # Computational homogenization
#
# ![](rve_homogenization.png)
#
# *Figure 1*: von Mises stress in an RVE with 5 stiff inclusions embedded in a softer matrix
# material that is loaded in shear. The problem is solved by using homogeneous Dirichlet
# boundary conditions (left) and (strong) periodic boundary conditions (right).
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`computational_homogenization.ipynb`](@__NBVIEWER_ROOT_URL__/examples/computational_homogenization.ipynb).
#-
#
# ## Introduction
#
# In this example we will solve the Representative Volume Element (RVE) problem for
# computational homogenization of linear elasticity and compute the effective/homogenized
# stiffness of an RVE with 5 stiff circular inclusions embedded in a softer matrix material
# (see Figure 1).
#
# It is possible to obtain upper and lower bounds on the stiffness analytically, see for
# example [Rule of mixtures](https://en.wikipedia.org/wiki/Rule_of_mixtures). An upper
# bound is obtained from the Voigt model, where the *strain* is assumed to be the same in
# the two constituents,
#
# ```math
# \mathsf{E}_\mathrm{Voigt} = v_\mathrm{m} \mathsf{E}_\mathrm{m} +
# (1 - v_\mathrm{m}) \mathsf{E}_\mathrm{i}
# ```
#
# where ``v_\mathrm{m}`` is the volume fraction of the matrix material, and where
# ``\mathsf{E}_\mathrm{m}`` and ``\mathsf{E}_\mathrm{i}`` are the individual stiffness for
# the matrix material and the inclusions, respectively. The lower bound is obtained from
# the Reuss model, where the *stress* is assumed to be the same in the two constituents,
#
# ```math
# \mathsf{E}_\mathrm{Reuss} = \left(v_\mathrm{m} \mathsf{E}_\mathrm{m}^{-1} +
# (1 - v_\mathrm{m}) \mathsf{E}_\mathrm{i}^{-1} \right)^{-1}.
# ```
#
# However, neither of these assumptions are, in general, very close to the "truth" which is
# why it is of interest to computationally find the homogenized properties for a given RVE.
#
# The canonical version of the RVE problem can be formulated as follows:
# For given homogenized field ``\bar{\boldsymbol{u}}``, ``\bar{\boldsymbol{\varepsilon}} =
# \boldsymbol{\varepsilon}[\bar{\boldsymbol{u}}]``, find ``\boldsymbol{u} \in
# \mathbb{U}_\Box``, ``\boldsymbol{t} \in \mathbb{T}_\Box`` such that
#
# ```math
# \frac{1}{|\Omega_\Box|} \int_{\Omega_\Box}\boldsymbol{\varepsilon}[\delta\boldsymbol{u}]
# : \mathsf{E} : \boldsymbol{\varepsilon}[\boldsymbol{u}]\ \mathrm{d}\Omega
# - \frac{1}{|\Omega_\Box|} \int_{\Gamma_\Box}\delta \boldsymbol{u} \cdot
# \boldsymbol{t}\ \mathrm{d}\Gamma = 0 \quad
# \forall \delta \boldsymbol{u} \in \mathbb{U}_\Box,\quad (1\mathrm{a})\\
# - \frac{1}{|\Omega_\Box|} \int_{\Gamma_\Box}\delta \boldsymbol{t} \cdot
# \boldsymbol{u}\ \mathrm{d}\Gamma = - \bar{\boldsymbol{\varepsilon}} :
# \left[ \frac{1}{|\Omega_\Box|} \int_{\Gamma_\Box}\delta \boldsymbol{t} \otimes
# [\boldsymbol{x} - \bar{\boldsymbol{x}}]\ \mathrm{d}\Gamma \right]
# \quad \forall \delta \boldsymbol{t} \in \mathbb{T}_\Box, \quad (1\mathrm{b})
# ```
#
# where ``\boldsymbol{u} = \bar{\boldsymbol{\varepsilon}} \cdot [\boldsymbol{x} -
# \bar{\boldsymbol{x}}] + \boldsymbol{u}^\mu``, where ``\Omega_\Box`` and ``|\Omega_\Box|``
# are the domain and volume of the RVE, where ``\Gamma_\Box`` is the boundary, and where
# ``\mathbb{U}_\Box``, ``\mathbb{T}_\Box`` are set of "sufficiently regular" functions
# defined on the RVE.
#
# This system is not solvable without introducing extra restrictions on ``\mathbb{U}_\Box``,
# ``\mathbb{T}_\Box``. In this example we will consider the common cases of Dirichlet
# boundary conditions and (strong) periodic boundary conditions.
#
# **Dirichlet boundary conditions**
#
# We can introduce the more restrictive sets of ``\mathbb{U}_\Box``:
#
# ```math
# \begin{align*}
# \mathbb{U}_\Box^\mathrm{D} &:= \left\{\boldsymbol{u} \in \mathbb{U}_\Box|\ \boldsymbol{u}
# = \bar{\boldsymbol{\varepsilon}} \cdot [\boldsymbol{x} - \bar{\boldsymbol{x}}]
# \ \mathrm{on}\ \Gamma_\Box\right\},\\
# \mathbb{U}_\Box^{\mathrm{D},0} &:= \left\{\boldsymbol{u} \in \mathbb{U}_\Box|\ \boldsymbol{u}
# = \boldsymbol{0}\ \mathrm{on}\ \Gamma_\Box\right\},
# \end{align*}
# ```
#
# and use these as trial and test sets to obtain a solvable RVE problem pertaining to
# Dirichlet boundary conditions. Eq. ``(1\mathrm{b})`` is trivially fulfilled, the second
# term of Eq. ``(1\mathrm{a})`` vanishes, and we are left with the following problem:
# Find ``\boldsymbol{u} \in \mathbb{U}_\Box^\mathrm{D}`` that solve
#
# ```math
# \frac{1}{|\Omega_\Box|} \int_{\Omega_\Box}\boldsymbol{\varepsilon}[\delta\boldsymbol{u}]
# : \mathsf{E} : \boldsymbol{\varepsilon}[\boldsymbol{u}]\ \mathrm{d}\Omega = 0
# \quad \forall \delta \boldsymbol{u} \in \mathbb{U}_\Box^{\mathrm{D},0}.
# ```
#
# Note that, since ``\boldsymbol{u} = \bar{\boldsymbol{\varepsilon}} \cdot [\boldsymbol{x} -
# \bar{\boldsymbol{x}}] + \boldsymbol{u}^\mu``, this problem is equivalent to solving for
# ``\boldsymbol{u}^\mu \in \mathbb{U}_\Box^{\mathrm{D},0}``, which is what we will do in
# the implementation.
#
# **Periodic boundary conditions**
#
# The RVE problem pertaining to periodic boundary conditions is obtained by restricting
# ``\boldsymbol{u}^\mu`` to be periodic, and ``\boldsymbol{t}`` anti-periodic across the
# RVE. Similarly as for Dirichlet boundary conditions, Eq. ``(1\mathrm{b})`` is directly
# fulfilled, and the second term in Eq. ``(1\mathrm{a})`` vanishes, with these restrictions,
# and we are left with the following problem:
# Find ``\boldsymbol{u}^\mu \in \mathbb{U}_\Box^{\mathrm{P},0}`` such that
#
# ```math
# \frac{1}{|\Omega_\Box|} \int_{\Omega_\Box}\boldsymbol{\varepsilon}[\delta\boldsymbol{u}]
# : \mathsf{E} : (\bar{\boldsymbol{\varepsilon}} + \boldsymbol{\varepsilon}
# [\boldsymbol{u}^\mu])\ \mathrm{d}\Omega = 0
# \quad \forall \delta \boldsymbol{u} \in \mathbb{U}_\Box^{\mathrm{P},0},
# ```
#
# where
#
# ```math
# \mathbb{U}_\Box^{\mathrm{P},0} := \left\{\boldsymbol{u} \in \mathbb{U}_\Box|
# \ \llbracket \boldsymbol{u} \rrbracket_\Box = \boldsymbol{0}
# \ \mathrm{on}\ \Gamma_\Box^+\right\}
# ```
#
# where ``\llbracket \bullet \rrbracket_\Box = \bullet(\boldsymbol{x}^+) -
# \bullet(\boldsymbol{x}^-)`` defines the "jump" over the RVE, i.e. the difference between
# the value on the image part ``\Gamma_\Box^+`` (coordinate ``\boldsymbol{x}^+``) and the
# mirror part ``\Gamma_\Box^-`` (coordinate ``\boldsymbol{x}^-``) of the boundary.
# To make sure this restriction holds in a strong sense we need a periodic mesh.
#
# Note that it would be possible to solve for the total ``\boldsymbol{u}`` directly by
# instead enforcing the jump to be equal to the jump in the macroscopic part,
# ``\boldsymbol{u}^\mathrm{M}``, i.e.
#
# ```math
# \llbracket \boldsymbol{u} \rrbracket_\Box =
# \llbracket \boldsymbol{u}^\mathrm{M} \rrbracket_\Box =
# \llbracket \bar{\boldsymbol{\varepsilon}} \cdot [\boldsymbol{x} - \bar{\boldsymbol{x}}]
# \rrbracket_\Box =
# \bar{\boldsymbol{\varepsilon}} \cdot [\boldsymbol{x}^+ - \boldsymbol{x}^-].
# ```
#
# **Homogenization of effective properties**
#
# In general it is necessary to compute the homogenized stress and the stiffness on the fly,
# but since we in this example consider linear elasticity it is possible to compute the
# effective properties once and for all for a given RVE configuration. We do this by
# computing sensitivity fields for every independent strain component (6 in 3D, 3 in 2D).
# Thus, for a 2D problem, as in the implementation below, we compute sensitivities
# ``\hat{\boldsymbol{u}}_{11}``, ``\hat{\boldsymbol{u}}_{22}``, and
# ``\hat{\boldsymbol{u}}_{12} = \hat{\boldsymbol{u}}_{21}`` by using
#
# ```math
# \bar{\boldsymbol{\varepsilon}} = \begin{pmatrix}1 & 0\\ 0 & 0\end{pmatrix}, \quad
# \bar{\boldsymbol{\varepsilon}} = \begin{pmatrix}0 & 0\\ 0 & 1\end{pmatrix}, \quad
# \bar{\boldsymbol{\varepsilon}} = \begin{pmatrix}0 & 0.5\\ 0.5 & 0\end{pmatrix}
# ```
#
# as the input to the RVE problem. When the sensitivies are solved we can compute the
# entries of the homogenized stiffness as follows
#
# ```math
# \mathsf{E}_{ijkl} = \frac{\partial\ \bar{\sigma}_{ij}}{\partial\ \bar{\varepsilon}_{kl}}
# = \bar{\sigma}_{ij}(\hat{\boldsymbol{u}}_{kl}),
# ```
#
# where the homogenized stress, ``\bar{\boldsymbol{\sigma}}(\boldsymbol{u})``, is computed
# as the volume average of the stress in the RVE, i.e.
#
# ```math
# \bar{\boldsymbol{\sigma}}(\boldsymbol{u}) :=
# \frac{1}{|\Omega_\Box|} \int_{\Omega_\Box} \boldsymbol{\sigma}\ \mathrm{d}\Omega =
# \frac{1}{|\Omega_\Box|} \int_{\Omega_\Box}
# \mathsf{E} : \boldsymbol{\varepsilon}[\boldsymbol{u}]\ \mathrm{d}\Omega.
# ```

# ## Commented program
#
# Now we will see how this can be implemented in `Ferrite`. What follows is a program
# with comments in between which describe the different steps.
#md # You can also find the same program without comments at the end of the page,
#md # see [Plain program](@ref homogenization-plain-program).

using Ferrite, SparseArrays, LinearAlgebra
using Test #src

# We first load the mesh file [`periodic-rve.msh`](periodic-rve.msh)
# ([`periodic-rve-coarse.msh`](periodic-rve-coarse.msh) for a coarser mesh). The mesh is
# generated with [`gmsh`](https://gmsh.info/), and we read it in as a `Ferrite` grid using
# the [`FerriteGmsh`](https://github.com/Ferrite-FEM/FerriteGmsh.jl) package:

using FerriteGmsh
#src notebook: use coarse mesh to decrease build time
#src   script: use the fine mesh
#src markdown: use the coarse mesh to decrease build time, but make it look like the fine
#nb ## grid = saved_file_to_grid("periodic-rve.msh")
#nb grid = saved_file_to_grid("periodic-rve-coarse.msh")
#jl ## grid = saved_file_to_grid("periodic-rve-coarse.msh")
#jl grid = saved_file_to_grid("periodic-rve.msh")
#md grid = saved_file_to_grid("periodic-rve.msh")
#-
#md grid = redirect_stdout(devnull) do                #hide
#md     saved_file_to_grid("periodic-rve-coarse.msh") #hide
#md end                                               #hide

grid = saved_file_to_grid("periodic-rve.msh") #src

# Next we construct the interpolation and quadrature rule, and combining them into
# cellvalues as usual:

dim = 2
ip = Lagrange{dim, RefTetrahedron, 1}()
qr = QuadratureRule{dim, RefTetrahedron}(2)
cellvalues = CellVectorValues(qr, ip);

# We define a dof handler with a displacement field `:u`:
dh = DofHandler(grid)
push!(dh, :u, 2)
close!(dh);

# Now we need to define boundary conditions. As discussed earlier we will solve the problem
# using (i) homogeneous Dirichlet boundary conditions, and (ii) periodic Dirichlet boundary
# conditions. We construct two different constraint handlers, one for each case. The
# [`Dirichlet`](@ref) boundary condition we have seen in many other examples. Here we simply
# define the condition that the field, `:u`, should have both components prescribed to `0`
# on the full boundary:

ch_dirichlet = ConstraintHandler(dh)
dirichlet = Dirichlet(
    :u,
    union(getfaceset.(Ref(grid), ["left", "right", "top", "bottom"])...),
    (x, t) ->  [0, 0],
    [1, 2]
)
add!(ch_dirichlet, dirichlet)
close!(ch_dirichlet)
update!(ch_dirichlet, 0.0)

# For periodic boundary conditions we use the [`PeriodicDirichlet`](@ref) constraint type,
# which is very similar to the `Dirichlet` type, but instead of a passing a faceset we pass
# a vector with "face pairs", i.e. the mapping between mirror and image parts of the
# boundary. In this example the `"left"` and `"bottom"` boundaries are mirrors, and the
# `"right"` and `"top"` boundaries are the mirrors.

ch_periodic = ConstraintHandler(dh);
periodic = PeriodicDirichlet(
    :u,
    ["left" => "right", "bottom" => "top"],
    [1, 2]
)
add!(ch_periodic, periodic)
close!(ch_periodic)
update!(ch_periodic, 0.0)

# This will now constrain any degrees of freedom located on the mirror boundaries to
# the matching degree of freedom on the image boundaries. Internally this will create
# a number of [`AffineConstraint`](@ref)s of the form `u_i = 1 * u_j + 0`:
# ```julia
# a = AffineConstraint(u_m, [u_i => 1], 0)
# ```
# where `u_m` is the degree of freedom on the mirror and `u_i` the matching one on the
# image part. `PeriodicDirichlet` is thus simply just a more convenient way of
# constructing such affine constraints since it computes the degree of freedom mapping
# automatically.
#
# To simplify things we group the constraint handlers into a named tuple

ch = (dirichlet = ch_dirichlet, periodic = ch_periodic);

# We can now construct the sparse matrix. Note that, since we are using affine constraints,
# which need to modify the matrix sparsity pattern in order to account for the constraint
# equations, we construct the matrix for the periodic case by passing both the dof handler
# and the constraint handler.

K = (
    dirichlet = create_sparsity_pattern(dh),
    periodic  = create_sparsity_pattern(dh, ch.periodic),
);

# We define the fourth order elasticity tensor for the matrix material, and define the
# inclusions to have 10 times higher stiffness

λ, μ = 1e10, 7e9 # Lamé parameters
δ(i,j) = i == j ? 1.0 : 0.0
Em = SymmetricTensor{4, 2}(
    (i,j,k,l) -> λ * δ(i,j) * δ(k,l) + μ * (δ(i,k) * δ(j,l) + δ(i,l) * δ(j,k))
)
Ei = 10 * Em;

# As mentioned above, in order to compute the apparent/homogenized stiffness we will solve
# the problem repeatedly with different macroscale strain tensors to compute the sensitvity
# of the homogenized stress, ``\bar{\boldsymbol{\sigma}}``, w.r.t. the macroscopic strain,
# ``\bar{\boldsymbol{\varepsilon}}``. The corresponding unit strains are defined below,
# and will result in three different right-hand-sides:

εᴹ = [
      SymmetricTensor{2,2}([1.0 0.0; 0.0 0.0]), # ε_11 loading
      SymmetricTensor{2,2}([0.0 0.0; 0.0 1.0]), # ε_22 loading
      SymmetricTensor{2,2}([0.0 0.5; 0.5 0.0]), # ε_12/ε_21 loading
];

# The assembly function is nothing strange, and in particular there is no impact from the
# choice of boundary conditions, so the same function can be used for both cases. Since
# we want to solve the system 3 times, once for each macroscopic strain component, we
# assemble 3 right-hand-sides.

function doassemble!(cellvalues::CellVectorValues, K::SparseMatrixCSC, dh::DofHandler, εᴹ)

    n_basefuncs = getnbasefunctions(cellvalues)
    ndpc = ndofs_per_cell(dh)
    Ke = zeros(ndpc, ndpc)
    fe = zeros(ndpc, length(εᴹ))
    f = zeros(ndofs(dh), length(εᴹ))
    assembler = start_assemble(K)

    for cell in CellIterator(dh)

        E = cellid(cell) in getcellset(dh.grid, "inclusions") ? Ei : Em
        reinit!(cellvalues, cell)
        fill!(Ke, 0)
        fill!(fe, 0)

        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            for i in 1:n_basefuncs
                δεi = shape_symmetric_gradient(cellvalues, q_point, i)
                for j in 1:n_basefuncs
                    δεj = shape_symmetric_gradient(cellvalues, q_point, j)
                    Ke[i, j] += (δεi ⊡ E ⊡ δεj) * dΩ
                end
                for (rhs, ε) in enumerate(εᴹ)
                    σᴹ = E ⊡ ε
                    fe[i, rhs] += ( - δεi ⊡ σᴹ) * dΩ
               end
            end
        end

        cdofs = celldofs(cell)
        assemble!(assembler, cdofs, Ke)
        f[cdofs, :] .+= fe
    end
    return f
end;

# We can now assemble the system. The assembly function modifies the matrix in-place, but
# return the right hand side(s) which we collect in another named tuple.

rhs = (
    dirichlet = doassemble!(cellvalues, K.dirichlet, dh, εᴹ),
    periodic  = doassemble!(cellvalues, K.periodic,  dh, εᴹ),
);

# The next step is to solve the systems. Since application of boundary conditions, using
# the [`apply!`](@ref) function, modifies both the matrix and the right hand sides we can
# not use it directly in this case since we want to reuse the matrix again for the next
# right hand sides. We could of course re-assemble the matrix for every right hand side,
# but that would not be very efficient. Instead we will use the [`get_rhs_data`](@ref)
# function, together with [`apply_rhs!`](@ref) in a later step. This will extract the
# necessary data from the matrix such that we can apply it for all the different right
# hand sides. Note that we call `apply!` with just the matrix and no right hand side.

rhsdata = (
    dirichlet = get_rhs_data(ch.dirichlet, K.dirichlet),
    periodic  = get_rhs_data(ch.periodic,  K.periodic),
)

apply!(K.dirichlet, ch.dirichlet)
apply!(K.periodic,  ch.periodic)

# We can now solve the problem(s). Note that we only use `apply_rhs!` in the loops below.
# The boundary conditions are already applied to the matrix above, so we only need to
# modify the right hand side.

u = (
    dirichlet = Vector{Float64}[],
    periodic  = Vector{Float64}[],
)

for i in 1:size(rhs.dirichlet, 2)
    rhs_i = @view rhs.dirichlet[:, i]                  # Extract this RHS
    apply_rhs!(rhsdata.dirichlet, rhs_i, ch.dirichlet) # Apply BC
    u_i = cholesky(Symmetric(K.dirichlet)) \ rhs_i     # Solve
    apply!(u_i, ch.dirichlet)                          # Apply BC on the solution
    push!(u.dirichlet, u_i)                            # Save the solution vector
end

for i in 1:size(rhs.periodic, 2)
    rhs_i = @view rhs.periodic[:, i]                   # Extract this RHS
    apply_rhs!(rhsdata.periodic, rhs_i, ch.periodic)   # Apply BC
    u_i = cholesky(Symmetric(K.periodic)) \ rhs_i      # Solve
    apply!(u_i, ch.periodic)                           # Apply BC on the solution
    push!(u.periodic, u_i)                             # Save the solution vector
end

# When the solution(s) are known we can compute the averaged stress,
# ``\bar{\boldsymbol{\sigma}}`` in the RVE. We define a function that does this, and also
# returns the von Mise stress in every quadrature point for visualization.

function compute_stress(cellvalues::CellVectorValues, dh::DofHandler, u, εᴹ)
    σvM_qpdata = zeros(getnquadpoints(cellvalues), getncells(dh.grid))
    σ̄Ω = zero(SymmetricTensor{2,2})
    Ω = 0.0 # Total volume
    for cell in CellIterator(dh)
        E = cellid(cell) in getcellset(dh.grid, "inclusions") ? Ei : Em
        reinit!(cellvalues, cell)
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            εμ = function_symmetric_gradient(cellvalues, q_point, u[celldofs(cell)])
            σ = E ⊡ (εᴹ + εμ)
            σvM_qpdata[q_point, cellid(cell)] = sqrt(3/2 * dev(σ) ⊡ dev(σ))
            Ω += dΩ # Update total volume
            σ̄Ω += σ * dΩ # Update integrated stress
        end
    end
    σ̄ = σ̄Ω / Ω
    return σvM_qpdata, σ̄
end;

# We now compute the homogenized stress and von Mise stress for all cases

σ̄ = (
    dirichlet = SymmetricTensor{2,2}[],
    periodic  = SymmetricTensor{2,2}[],
)
σ = (
     dirichlet = Vector{Float64}[],
     periodic  = Vector{Float64}[],
)

projector = L2Projector(ip, grid)

for i in 1:3
    σ_qp, σ̄_i = compute_stress(cellvalues, dh, u.dirichlet[i], εᴹ[i])
    proj = project(projector, σ_qp, qr; project_to_nodes=false)
    push!(σ.dirichlet, proj)
    push!(σ̄.dirichlet, σ̄_i)
end

for i in 1:3
    σ_qp, σ̄_i = compute_stress(cellvalues, dh, u.periodic[i], εᴹ[i])
    proj = project(projector, σ_qp, qr; project_to_nodes=false)
    push!(σ.periodic, proj)
    push!(σ̄.periodic, σ̄_i)
end

# The remaining thing is to compute the homogenized stiffness. As mentioned in the
# introduction we can find all the components from the average stress of the sensitivity
# fields that we have solved for
#
# ```math
# \mathsf{E}_{ijkl} = \bar{\sigma}_{ij}(\hat{\boldsymbol{u}}_{kl}).
# ```
#
# So we have now already computed all the components, and just need to gather the data in
# a fourth order tensor:

E_dirichlet = SymmetricTensor{4,2}((i, j, k, l) -> begin
    if k == l == 1
        σ̄.dirichlet[1][i, j] # ∂σ∂ε_**11
    elseif k == l == 2
        σ̄.dirichlet[2][i, j] # ∂σ∂ε_**22
    else
        σ̄.dirichlet[3][i, j] # ∂σ∂ε_**12 and ∂σ∂ε_**21
    end
end)

E_periodic = SymmetricTensor{4,2}((i, j, k, l) -> begin
    if k == l == 1
        σ̄.periodic[1][i, j]
    elseif k == l == 2
        σ̄.periodic[2][i, j]
    else
        σ̄.periodic[3][i, j]
    end
end);

# We can check that the result are what we expect, namely that the stiffness with Dirichlet
# boundary conditions is higher than when using periodic boundary conditions, and that
# the Reuss assumption is an lower bound, and the Voigt assumption a upper bound. We first
# compute the volume fraction of the matrix, and then the Voigt and Reuss bounds:

function matrix_volume_fraction(grid, cellvalues)
    V  = 0.0 # Total volume
    Vm = 0.0 # Volume of the matrix
    for c in CellIterator(grid)
        reinit!(cellvalues, c)
        is_matrix = !(cellid(c) in getcellset(grid, "inclusions"))
        for qp in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, qp)
            V += dΩ
            if is_matrix
                Vm += dΩ
            end
        end
    end
    return Vm / V
end

vm = matrix_volume_fraction(grid, cellvalues)
#-
E_voigt = vm * Em + (1-vm) * Ei
E_reuss = inv(vm * inv(Em) + (1-vm) * inv(Ei));

# We can now compare the different computed stiffness tensors. We expect
# ``E_\mathrm{Reuss} \leq E_\mathrm{PeriodicBC} \leq E_\mathrm{DirichletBC} \leq
# E_\mathrm{Voigt}``. A simple thing to compare are the eigenvalues of the tensors. Here
# we look at the first eigenvalue:

ev = (first ∘ eigvals).((E_reuss, E_periodic, E_dirichlet, E_voigt))
@test issorted(ev) #src
round.(ev; digits=-8)

# Finally, we export the solution and the stress field to a VTK file. For the export we
# also compute the macroscopic part of the displacement.

chM = ConstraintHandler(dh)
add!(chM, Dirichlet(:u, Set(1:getnnodes(grid)), (x, t) -> εᴹ[Int(t)] ⋅ x, [1, 2]))
close!(chM)
uM = zeros(ndofs(dh))

vtk_grid("homogenization", dh) do vtk
    for i in 1:3
        ## Compute macroscopic solution
        update!(chM, i)
        apply!(uM, chM)
        ## Dirichlet
        vtk_point_data(vtk, dh, uM + u.dirichlet[i], "_dirichlet_$i")
        vtk_point_data(vtk, projector, σ.dirichlet[i], "σvM_dirichlet_$i")
        ## Periodic
        vtk_point_data(vtk, dh, uM + u.periodic[i], "_periodic_$i")
        vtk_point_data(vtk, projector, σ.periodic[i], "σvM_periodic_$i")
    end
end;

# Just another way to compute the stiffness for testing purposes               #src
function homogenize_test(u::Matrix, dh, cv, E_incl, E_mat)                     #src
    ĒΩ = zero(SymmetricTensor{4,2})                                            #src
    Ω = 0.0                                                                    #src
    ue = zeros(ndofs_per_cell(dh), 3)                                          #src
    for cell in CellIterator(dh)                                               #src
        reinit!(cv, cell)                                                      #src
        for (localdof, globaldof) in enumerate(celldofs(cell))                 #src
            for i in 1:3                                                       #src
                ue[localdof, i] = u[globaldof, i]                              #src
            end                                                                #src
        end                                                                    #src
        E = cellid(cell) in getcellset(dh.grid, "inclusions") ? E_incl : E_mat #src
        for qp in 1:getnquadpoints(cv)                                         #src
            dΩ = getdetJdV(cv, qp)                                             #src
            Ω += dΩ                                                            #src
            ## compute u^ij and u^kl                                           #src
            Ē′ = SymmetricTensor{4,2}((i, j, k, l) -> begin                    #src
                ij = i == j == 1 ? 1 : i == j == 2 ? 2 : 3                     #src
                kl = k == l == 1 ? 1 : k == l == 2 ? 2 : 3                     #src
                εij = function_symmetric_gradient(cv, qp, view(ue, :, ij)) +   #src
                          symmetric((basevec(Vec{2}, i) ⊗ basevec(Vec{2}, j))) #src
                εkl = function_symmetric_gradient(cv, qp, view(ue, :, kl)) +   #src
                          symmetric((basevec(Vec{2}, k) ⊗ basevec(Vec{2}, l))) #src
                return (εij ⊡ E ⊡ εkl) * dΩ                                    #src
            end)                                                               #src
            ĒΩ += Ē′                                                           #src
        end                                                                    #src
    end                                                                        #src
    return ĒΩ / Ω                                                              #src
end                                                                            #src

@test homogenize_test(reduce(hcat, u.dirichlet), dh, cellvalues, Ei, Em) ≈ E_dirichlet #src
@test homogenize_test(reduce(hcat, u.periodic), dh, cellvalues, Ei, Em) ≈ E_periodic #src

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


text:
# Version management

export version

version() = Base.libllvm_version

function runtime_version()
    # FIXME: add a proper C API to LLVM
    version_print = unsafe_string(
        ccall((:_ZN4llvm16LTOCodeGenerator16getVersionStringEv, libllvm), Cstring, ()))
    m = match(r"LLVM version (?<version>.+)", version_print)
    m === nothing && error("Unrecognized version string: '$version_print'")
    if endswith(m[:version], "jl")
        # strip the "jl" SONAME suffix (JuliaLang/julia#33058)
        # (LLVM does never report a prerelease version anyway)
        VersionNumber(m[:version][1:end-2])
    else
        VersionNumber(m[:version])
    end
end

function is_asserts()
    Libdl.dlopen(libllvm) do handle
        Libdl.dlsym(handle, "_ZN4llvm23EnableABIBreakingChecksE"; throw_error=false) !== nothing
    end
end

ext: jl


text:
""" Free diagrams in a category.

A [free diagram](https://ncatlab.org/nlab/show/free+diagram) in a category is a
diagram whose shape is a free category. Examples include the empty diagram,
pairs of objects, discrete diagrams, parallel morphisms, spans, and cospans.
Limits and colimits are most commonly taken over free diagrams.
"""
module FreeDiagrams
export AbstractFreeDiagram, FreeDiagram, FixedShapeFreeDiagram, DiscreteDiagram,
  EmptyDiagram, ObjectPair, Span, Cospan, Multispan, Multicospan,
  SMultispan, SMulticospan, ParallelPair, ParallelMorphisms,
  ob, hom, dom, codom, apex, legs, feet, left, right,
  nv, ne, src, tgt, vertices, edges, has_vertex, has_edge,
  add_vertex!, add_vertices!, add_edge!, add_edges!,
  SquareDiagram, top, bottom

using AutoHashEquals
using StaticArrays: StaticVector, SVector, @SVector

using ...Present, ...Theories, ...CSetDataStructures, ...Graphs
import ...Theories: ob, hom, dom, codom, top, bottom, left, right
using ...Graphs.BasicGraphs: TheoryGraph

# Diagrams of fixed shape
#########################

""" Abstract type for free diagram of fixed shape.
"""
abstract type FixedShapeFreeDiagram{Ob} end

""" Discrete diagram: a diagram whose only morphisms are identities.
"""
@auto_hash_equals struct DiscreteDiagram{Ob,Objects<:AbstractVector{Ob}} <:
    FixedShapeFreeDiagram{Ob}
  objects::Objects
end

const EmptyDiagram{Ob} = DiscreteDiagram{Ob,<:StaticVector{0}}
const ObjectPair{Ob} = DiscreteDiagram{Ob,<:StaticVector{2}}

EmptyDiagram{Ob}() where Ob = DiscreteDiagram(@SVector Ob[])
ObjectPair(first, second) = DiscreteDiagram(SVector(first, second))

ob(d::DiscreteDiagram) = d.objects

Base.iterate(d::DiscreteDiagram, args...) = iterate(d.objects, args...)
Base.eltype(d::DiscreteDiagram) = eltype(d.objects)
Base.length(d::DiscreteDiagram) = length(d.objects)
Base.getindex(d::DiscreteDiagram, i) = d.objects[i]
Base.firstindex(d::DiscreteDiagram) = firstindex(d.objects)
Base.lastindex(d::DiscreteDiagram) = lastindex(d.objects)

""" Multispan of morphisms in a category.

A [multispan](https://ncatlab.org/nlab/show/multispan) is like a [`Span`](@ref)
except that it may have a number of legs different than two. A colimit of this
shape is a pushout.
"""
@auto_hash_equals struct Multispan{Ob,Legs<:AbstractVector} <:
    FixedShapeFreeDiagram{Ob}
  apex::Ob
  legs::Legs
end

function Multispan(legs::AbstractVector)
  !isempty(legs) || error("Empty list of legs but no apex given")
  allequal(dom.(legs)) || error("Legs $legs do not have common domain")
  Multispan(dom(first(legs)), legs)
end

const SMultispan{N,Ob} = Multispan{Ob,<:StaticVector{N}}

SMultispan{N}(apex, legs::Vararg{Any,N}) where N =
  Multispan(apex, SVector{N}(legs...))
SMultispan{0}(apex) = Multispan(apex, SVector{0,typeof(id(apex))}())
SMultispan{N}(legs::Vararg{Any,N}) where N = Multispan(SVector{N}(legs...))

""" Span of morphims in a category.

A common special case of [`Multispan`](@ref). See also [`Cospan`](@ref).
"""
const Span{Ob} = SMultispan{2,Ob}

apex(span::Multispan) = span.apex
legs(span::Multispan) = span.legs
feet(span::Multispan) = map(codom, span.legs)
left(span::Span) = span.legs[1]
right(span::Span) = span.legs[2]

Base.iterate(span::Multispan, args...) = iterate(span.legs, args...)
Base.eltype(span::Multispan) = eltype(span.legs)
Base.length(span::Multispan) = length(span.legs)

""" Multicospan of morphisms in a category.

A multicospan is like a [`Cospan`](@ref) except that it may have a number of
legs different than two. A limit of this shape is a pullback.
"""
@auto_hash_equals struct Multicospan{Ob,Legs<:AbstractVector} <:
    FixedShapeFreeDiagram{Ob}
  apex::Ob
  legs::Legs
end

function Multicospan(legs::AbstractVector)
  !isempty(legs) || error("Empty list of legs but no apex given")
  allequal(codom.(legs)) || error("Legs $legs do not have common codomain")
  Multicospan(codom(first(legs)), legs)
end

const SMulticospan{N,Ob} = Multicospan{Ob,<:StaticVector{N}}

SMulticospan{N}(apex, legs::Vararg{Any,N}) where N =
  Multicospan(apex, SVector{N}(legs...))
SMulticospan{0}(apex) = Multicospan(apex, SVector{0,typeof(id(apex))}())
SMulticospan{N}(legs::Vararg{Any,N}) where N = Multicospan(SVector{N}(legs...))

""" Cospan of morphisms in a category.

A common special case of [`Multicospan`](@ref). See also [`Span`](@ref).
"""
const Cospan{Ob} = SMulticospan{2,Ob}

apex(cospan::Multicospan) = cospan.apex
legs(cospan::Multicospan) = cospan.legs
feet(cospan::Multicospan) = map(dom, cospan.legs)
left(cospan::Cospan) = cospan.legs[1]
right(cospan::Cospan) = cospan.legs[2]

Base.iterate(cospan::Multicospan, args...) = iterate(cospan.legs, args...)
Base.eltype(cospan::Multicospan) = eltype(cospan.legs)
Base.length(cospan::Multicospan) = length(cospan.legs)

""" Parallel morphims in a category.

[Parallel morphisms](https://ncatlab.org/nlab/show/parallel+morphisms) are just
morphisms with the same domain and codomain. A (co)limit of this shape is a
(co)equalizer.

For the common special case of two morphisms, see [`ParallelPair`](@ref).
"""
@auto_hash_equals struct ParallelMorphisms{Ob,Homs<:AbstractVector} <:
    FixedShapeFreeDiagram{Ob}
  dom::Ob
  codom::Ob
  homs::Homs
end

function ParallelMorphisms(homs::AbstractVector)
  @assert !isempty(homs) && allequal(dom.(homs)) && allequal(codom.(homs))
  ParallelMorphisms(dom(first(homs)), codom(first(homs)), homs)
end

""" Pair of parallel morphisms in a category.

A common special case of [`ParallelMorphisms`](@ref).
"""
const ParallelPair{Ob} = ParallelMorphisms{Ob,<:StaticVector{2}}

function ParallelPair(first, last)
  dom(first) == dom(last) ||
    error("Domains of parallel pair do not match: $first vs $last")
  codom(first) == codom(last) ||
    error("Codomains of parallel pair do not match: $first vs $last")
  ParallelMorphisms(dom(first), codom(first), SVector(first, last))
end

dom(para::ParallelMorphisms) = para.dom
codom(para::ParallelMorphisms) = para.codom
hom(para::ParallelMorphisms) = para.homs

Base.iterate(para::ParallelMorphisms, args...) = iterate(para.homs, args...)
Base.eltype(para::ParallelMorphisms) = eltype(para.homs)
Base.length(para::ParallelMorphisms) = length(para.homs)
Base.getindex(para::ParallelMorphisms, i) = para.homs[i]
Base.firstindex(para::ParallelMorphisms) = firstindex(para.homs)
Base.lastindex(para::ParallelMorphisms) = lastindex(para.homs)

allequal(xs::AbstractVector) = all(isequal(x, xs[1]) for x in xs[2:end])

# Commutative squares
#####################

"""    SquareDiagram(top, bottom, left, right)

creates a square diagram in a category, which forms the 2-cells of the double category Sq(C).
The four 1-cells are given in top, bottom, left, right order, to match the GAT of a double category.
"""
@auto_hash_equals struct SquareDiagram{Ob,Homs<:AbstractVector} <:
    FixedShapeFreeDiagram{Ob}
  corners::Vector{Ob}
  sides::Homs
end

function SquareDiagram(homs::AbstractVector)
  length(homs) == 4 || error("Square diagrams accept exactly 4 homs, in order top, bottom, left, right")
  obs = [dom(homs[1]), dom(homs[2]), dom(homs[4]), codom(homs[4])]

  # checking well-formedness
  # top-left share domains
  @assert dom(homs[1]) == dom(homs[3])
  # bottom-right share codomains
  @assert codom(homs[2]) == codom(homs[4])
  # left-bottom intersection
  @assert codom(homs[3]) == dom(homs[2])
  # top-right intersection
  @assert codom(homs[1]) == dom(homs[4])

  SquareDiagram(obs, homs)
end

SquareDiagram(top, bottom, left, right) = SquareDiagram([top, bottom, left, right])

ob(sq::SquareDiagram) = sq.corners
hom(sq::SquareDiagram) = sq.sides

top(sq::SquareDiagram) = sq.sides[1]
bottom(sq::SquareDiagram) = sq.sides[2]
left(sq::SquareDiagram) = sq.sides[3]
right(sq::SquareDiagram) = sq.sides[4]

# General diagrams
##################

@present TheoryFreeDiagram <: TheoryGraph begin
  Ob::Data
  Hom::Data
  ob::Attr(V,Ob)
  hom::Attr(E,Hom)
end

const FreeDiagram = ACSetType(TheoryFreeDiagram, index=[:src,:tgt])

# XXX: This is needed because we cannot control the supertype of C-set types.
const _AbstractFreeDiagram = AbstractACSetType(TheoryFreeDiagram)
const AbstractFreeDiagram{Ob} =
  Union{_AbstractFreeDiagram{Ob},FixedShapeFreeDiagram{Ob}}

ob(d::FreeDiagram, args...) = subpart(d, args..., :ob)
hom(d::FreeDiagram, args...) = subpart(d, args..., :hom)

function FreeDiagram(obs::Vector{Ob},
                     homs::Vector{Tuple{Int,Int,Hom}}) where {Ob,Hom}
  @assert all(obs[s] == dom(f) && obs[t] == codom(f) for (s,t,f) in homs)
  d = FreeDiagram{Ob,Hom}()
  add_vertices!(d, length(obs), ob=obs)
  add_edges!(d, getindex.(homs,1), getindex.(homs,2), hom=last.(homs))
  return d
end

# Conversion of fixed shapes
#---------------------------

function FreeDiagram(discrete::DiscreteDiagram{Ob}) where Ob
  d = FreeDiagram{Ob,Nothing}()
  add_vertices!(d, length(discrete), ob=collect(discrete))
  return d
end

function FreeDiagram(span::Multispan{Ob}) where Ob
  d = FreeDiagram{Ob,eltype(span)}()
  v0 = add_vertex!(d, ob=apex(span))
  vs = add_vertices!(d, length(span), ob=feet(span))
  add_edges!(d, fill(v0, length(span)), vs, hom=legs(span))
  return d
end

function FreeDiagram(cospan::Multicospan{Ob}) where Ob
  d = FreeDiagram{Ob,eltype(cospan)}()
  vs = add_vertices!(d, length(cospan), ob=feet(cospan))
  v0 = add_vertex!(d, ob=apex(cospan))
  add_edges!(d, vs, fill(v0, length(cospan)), hom=legs(cospan))
  return d
end

function FreeDiagram(para::ParallelMorphisms{Ob}) where Ob
  d = FreeDiagram{Ob,eltype(para)}()
  add_vertices!(d, 2, ob=[dom(para), codom(para)])
  add_edges!(d, fill(1,length(para)), fill(2,length(para)), hom=hom(para))
  return d
end

function FreeDiagram(sq::SquareDiagram{Ob}) where Ob
    top, bottom, left, right = sq.sides
    # check that the domains and codomains match
    #   1   -top->   3
    #   |            |
    # left         right
    #   v            v
    #   2  -bottom-> 4

    @assert codom(top) == dom(right)
    @assert dom(top) == dom(left)
    @assert dom(bottom) == codom(left)
    @assert codom(bottom) == codom(right)

    V = [dom(left), codom(left), dom(right), codom(right)]
    E = [(1,3, top), (2,4, bottom), (1,2, left), (3,4, right)]
    return FreeDiagram(V, E)
end
end


text:
"""
    rescale_hidden!(rbm, λ::AbstractArray)

For continuous hidden units with a scale parameter,
scales parameters such that hidden unit activations
are divided by `λ`. For other hidden units does
nothing. The resulting RBM is equivalent to the
original one.
"""
function rescale_hidden!(rbm::RBM, λ::AbstractArray)
    @assert size(rbm.hidden) == size(λ)
    if rescale_activations!(rbm.hidden, λ)
        rbm.w .*= reshape(λ, map(one, size(rbm.visible))..., size(rbm.hidden)...)
    end
    return rbm
end

"""
    rescale_activations!(layer, λ::AbstractArray)

For continuous layers with scale parameters, re-parameterizes
such that unit activations are divided by `λ`, and returns `true`.
For other layers just returns `false`.
"""
rescale_activations!(layer::Union{Binary,Spin,Potts}, λ::AbstractArray) = false

#= Note that λ < 0 can lead to trouble, e.g. for ReLU which
must have positive activations. So we dissallow it below. =#

function rescale_activations!(layer::Union{Gaussian,ReLU}, λ::AbstractArray)
    @assert size(layer) == size(λ)
    @assert all(λ .> 0)
    layer.θ .*= λ
    layer.γ .*= λ.^2
    return true
end

function rescale_activations!(layer::dReLU, λ::AbstractArray)
    @assert size(layer) == size(λ)
    @assert all(λ .> 0)
    layer.θp .*= λ
    layer.θn .*= λ
    layer.γp .*= λ.^2
    layer.γn .*= λ.^2
    return true
end

function rescale_activations!(layer::Union{pReLU,xReLU}, λ::AbstractArray)
    @assert size(layer) == size(λ)
    @assert all(λ .> 0) # makes life simpler
    layer.θ .*= λ
    layer.Δ .*= λ
    layer.γ .*= λ.^2
    return true
end


text:
# This file is a part of Julia, but is derived from
# https://github.com/floitsch/double-conversion which has the following license
#
# Copyright 2006-2014, the V8 project authors. All rights reserved.
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
#     * Redistributions of source code must retain the above copyright
#       notice, this list of conditions and the following disclaimer.
#     * Redistributions in binary form must reproduce the above
#       copyright notice, this list of conditions and the following
#       disclaimer in the documentation and/or other materials provided
#       with the distribution.
#     * Neither the name of Google Inc. nor the names of its
#       contributors may be used to endorse or promote products derived
#       from this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
# A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
# OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

function normalizedexponent(significand, exponent::Int32)
    significand = UInt64(significand)
    while (significand & HiddenBit(Float64)) == 0
        significand <<= UInt64(1)
        exponent -= Int32(1)
    end
    return exponent
end

function bignumdtoa(v,mode,requested_digits::Int,buffer,bignums)
    significand = _significand(v)
    exponent = _exponent(v)
    lower_boundary_is_closer = lowerboundaryiscloser(v)
    need_boundary_deltas = mode == SHORTEST

    is_even = (significand & 1) == 0
    normalized_exponent = normalizedexponent(significand, exponent)
    estimated_power = estimatepower(Int(normalized_exponent))

    if mode == FIXED && -estimated_power - 1 > requested_digits
        buffer[1] = 0
        len = 1
        decimal_point = -requested_digits
        return true, len, decimal_point
    end
    num, den, minus, plus = bignums[1], bignums[2], bignums[3], bignums[4]
    initialscaledstartvalues!(significand,exponent,lower_boundary_is_closer,
                              estimated_power,need_boundary_deltas,
                              num,den,minus,plus)
    decimal_point = fixupmultiply10!(estimated_power,is_even,num,den,minus,plus)
    if mode == SHORTEST
        len = generateshortestdigits!(num,den,minus,plus,is_even,buffer)
    elseif mode == FIXED
        len, decimal_point = bignumtofixed!(requested_digits,num,den,buffer,decimal_point)
    elseif mode == PRECISION
        len, decimal_point = generatecounteddigits!(requested_digits,num,den,buffer,decimal_point)
    end
    buffer[len] = 0
    return true, len, decimal_point
end

function generateshortestdigits!(num,den,minus,plus,is_even,buffer)
    minus == plus && (plus = minus)
    len = 1
    while true
        digit = Bignums.dividemodulointbignum!(num,den)
        buffer[len] = 0x30 + (digit % UInt8)
        len += 1
        in_delta_room_minus = is_even ?
            Bignums.lessequal(num,minus) : Bignums.less(num,minus)
        in_delta_room_plus = is_even ?
            Bignums.pluscompare(num,plus,den) >= 0: Bignums.pluscompare(num,plus,den) > 0
        if !in_delta_room_minus && !in_delta_room_plus
            Bignums.times10!(num)
            Bignums.times10!(minus)
            minus != plus && Bignums.times10!(plus)
        elseif in_delta_room_minus && in_delta_room_plus
            compare = Bignums.pluscompare(num,num,den)
            if compare < 0
            elseif compare > 0
                buffer[len - 1] += 1
            else
                if (buffer[len - 1] - 0x30) % 2 == 0
                else
                    buffer[len - 1] += 1
                end
            end
            return len
        elseif in_delta_room_minus
            return len
        else
            buffer[len - 1] += 1
            return len
        end
    end
end

function generatecounteddigits!(count,num,den,buffer,decimal_point)
    for i = 1:(count-1)
        digit = Bignums.dividemodulointbignum!(num,den)
        buffer[i] = 0x30 + (digit % UInt8)
        Bignums.times10!(num)
    end
    digit = Bignums.dividemodulointbignum!(num,den)
    if Bignums.pluscompare(num,num,den) >= 0
        digit += 1
    end
    buffer[count] = 0x30 + (digit % UInt8)
    for i = count:-1:2
        buffer[i] != 0x30 + 10 && break
        buffer[i] = 0x30
        buffer[i - 1] += 1
    end
    if buffer[1] == 0x30 + 10
        buffer[1] = 0x31
        decimal_point += 1
    end
    len = count+1
    return len, decimal_point
end

function bignumtofixed!(requested_digits,num,den,buffer,decimal_point)
    if -decimal_point > requested_digits
        decimal_point = -requested_digits
        len = 1
        return len, decimal_point
    elseif -decimal_point == requested_digits
        Bignums.times10!(den)
        if Bignums.pluscompare(num,num,den) >= 0
            buffer[1] = 0x31
            len = 2
            decimal_point += 1
        else
            len = 1
        end
        return len, decimal_point
    else
        needed_digits = decimal_point + requested_digits
        len, decimal_point = generatecounteddigits!(
              needed_digits,num,den,buffer,decimal_point)
    end
    return len, decimal_point
end

const k1Log10 = 0.30102999566398114
const kSignificandSize = SignificandSize(Float64)
estimatepower(exponent::Int) = ceil(Int,(exponent + kSignificandSize - 1) * k1Log10 - 1e-10)

function init3!(
        significand,exponent,estimated_power,need_boundary_deltas,
        num,den,minus,plus)
    Bignums.assignuint64!(num,UInt64(significand))
    Bignums.shiftleft!(num,exponent)
    Bignums.assignpoweruint16!(den,UInt16(10),estimated_power)
    if need_boundary_deltas
        Bignums.shiftleft!(den,1)
        Bignums.shiftleft!(num,1)
        Bignums.assignuint16!(plus,UInt16(1))
        Bignums.shiftleft!(plus,exponent)
        Bignums.assignuint16!(minus,UInt16(1))
        Bignums.shiftleft!(minus,exponent)
    else
        Bignums.zero!(plus)
        Bignums.zero!(minus)
    end
    return
end

function init1!(
        significand,exponent,estimated_power,need_boundary_deltas,
        num,den,minus,plus)
    Bignums.assignuint64!(num,UInt64(significand))
    Bignums.assignpoweruint16!(den,UInt16(10),estimated_power)
    Bignums.shiftleft!(den,-exponent)
    if need_boundary_deltas
        Bignums.shiftleft!(den,1)
        Bignums.shiftleft!(num,1)
        Bignums.assignuint16!(plus,UInt16(1))
        Bignums.assignuint16!(minus,UInt16(1))
    else
        Bignums.zero!(plus)
        Bignums.zero!(minus)
    end
    return
end

function init2!(
        significand,exponent,estimated_power,need_boundary_deltas,
        num,den,minus,plus)
    power_ten = num
    Bignums.assignpoweruint16!(power_ten,UInt16(10),-estimated_power)
    if need_boundary_deltas
        Bignums.assignbignum!(plus,power_ten)
        Bignums.assignbignum!(minus,power_ten)
    else
        Bignums.zero!(plus)
        Bignums.zero!(minus)
    end
    Bignums.multiplybyuint64!(num,UInt64(significand))
    Bignums.assignuint16!(den,UInt16(1))
    Bignums.shiftleft!(den,-exponent)
    if need_boundary_deltas
        Bignums.shiftleft!(num,1)
        Bignums.shiftleft!(den,1)
    end
    return
end

function initialscaledstartvalues!(significand,
            exponent,lower_boundary_is_closer,estimated_power,
            need_boundary_deltas,num,den,minus,plus)
    if exponent >= 0
        init3!(significand, exponent, estimated_power, need_boundary_deltas,num,den,minus,plus)
    elseif estimated_power >= 0
        init1!(significand, exponent, estimated_power, need_boundary_deltas,num,den,minus,plus)
    else
        init2!(significand, exponent, estimated_power, need_boundary_deltas,num,den,minus,plus)
    end
    if need_boundary_deltas && lower_boundary_is_closer
        Bignums.shiftleft!(den,1)
        Bignums.shiftleft!(num,1)
        Bignums.shiftleft!(plus,1)
    end
    return
end

function fixupmultiply10!(estimated_power,is_even,num,den,minus,plus)
    in_range = is_even ? Bignums.pluscompare(num,plus,den) >= 0 :
                         Bignums.pluscompare(num,plus,den) > 0
    if in_range
        decimal_point = estimated_power + 1
    else
        decimal_point = estimated_power
        Bignums.times10!(num)
        if minus == plus
            Bignums.times10!(minus)
            Bignums.assignbignum!(plus,minus)
        else
            Bignums.times10!(minus)
            Bignums.times10!(plus)
        end
    end
    return decimal_point
end


text:
include("layout.jl")                    # layout
include(joinpath("api", "endpoint.jl")) # Endpoint
include("wasm.jl")                      # WASM
include("front.jl")                     # Front

if PROGRAM_FILE == basename(@__FILE__)

using Bukdu # pipeline Conn routes resources get plug Router
import .Bukdu.Actions: index, show, new, edit, create, delete, update
using .Front # WelcomeController
using .Endpoint # CustomerController
using .WASM: WasmController, hello_js, hello_wast
using Sockets

pipeline(:api) do conn::Conn
end

routes(:api) do
    resources("/customers", CustomerController, :id=>Int)
end

routes(:wasm) do
    get("/wasm", WasmController, index)
    get("/hello.js", WasmController, hello_js)
    get("/hello.wast", WasmController, hello_wast)
    plug(Plug.Static, at="/", from=normpath(@__DIR__, "public"))
end

routes(:front) do
    get("/", WelcomeController, index)
end

Bukdu.start(Front.server.port; host=Front.server.host)

Router.call(get, "/") #
# CLI.routes()

Base.JLOptions().isinteractive==0 && wait()

# Bukdu.stop()

end # if


text:
"""
`init_subspec!(m::AnSchorfheide)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::AnSchorfheide)
    if subspec(m) == "ss0"
        return
    else
        error("This subspec should be a 0")
    end
end


text:
module test_loltools_matchv5

using Test
using LOLTools.MatchV5

@test MatchV5.matches_by_puuid isa Function

include("helpers.jl")

struct MatchController <: ApplicationController
    conn::Conn
end

dict_matches_ids = Dict{String, Vector{String}}(
    "puuid1" => String[],
    "puuid2" => String["KR_12345", "KR_12346"],
)

function matches_by_puuid(c::MatchController)
    puuid = c.params.puuid
    matches_ids = dict_matches_ids[puuid]
    render(JSON, matches_ids)
end

res_match = (metadata=0, info=1)

function match_by_match_id(c::MatchController)
    render(JSON, res_match)
end

res_match_timeline = (metadata=0, info=1)

function match_timeline_by_match_id(c::MatchController)
    render(JSON, res_match_timeline)
end

routes() do
    get("/lol/match/v5/matches/by-puuid/:puuid/ids", MatchController, matches_by_puuid)
	get("/lol/match/v5/matches/:matchId", MatchController, match_by_match_id)
	get("/lol/match/v5/matches/:matchId/timeline", MatchController, match_timeline_by_match_id)
end

api_key = ""
region = "asia"

puuid = "puuid1"
@test MatchV5.matches_by_puuid(api_key, region, puuid; action=mock_action) == String[]

puuid = "puuid2"
@test MatchV5.matches_by_puuid(api_key, region, puuid; action=mock_action) == String["KR_12345", "KR_12346"]

matchId = "KR_12345"
MatchV5.match_by_match_id(api_key, region, matchId; action=mock_action)
MatchV5.match_timeline_by_match_id(api_key, region, matchId; action=mock_action)

end # module test_loltools_matchv5


text:
@testset "Summarizers" begin
    using Test

    # requires Base NamedTuples and DataStreams/jq/gangy

    using DataStreams
    using StreamModels

    using StreamModels: update!, summarize, CategoricalSummarizer, make_summarizer

    source = Data.Table((a = collect(1:10),
                         b = rand(10),
                         c = repeat(["a", "b"], inner=5)))
    sch = Data.schema(source)
    types = Dict(Symbol(k) => Data.types(sch)[sch[k]] for k in Data.header(sch))

    @testset "CategoricalSummarizer computes unique values" begin
        iter = Data.rows(source)

        summ = StreamModels.make_summarizer(StreamModels.Terms.Eval(:c), types)

        for nt in iter
            update!(summ, nt)
        end

        @test unique(source[:c]) == summ.uniq
    end

    @testset "If no summaries are needed, none are calculated" begin
        sch = summarize(source, @formula(~ 1 + a * b))
        @test !haskey(sch.metadata, :summaries)
    end

    @testset "One summarizer per unique variable" begin
        f = @formula( ~ 1 + a + b + c + c)
        summarizers = unique(reduce(vcat, [], make_summarizer(t, types) for t in f.term.rhs))
        @test length(summarizers) == 1

        f = @formula( ~ 1 + a*b*c)
        summarizers = unique(reduce(vcat, [], make_summarizer(t, types) for t in f.term.rhs))
        @test length(summarizers) == 1
    end

    @testset "Unique values stored in sch.metadata for categorical variables" begin
        sch = summarize(source, @formula( ~ c))
        @test collect(keys(sch.metadata[:summaries])) == Any[:c]
        @test collect(keys(sch.metadata[:summaries][:c])) == Any[:unique]
        @test sch.metadata[:summaries][:c][:unique] == unique(source[:c])
    end
end #testset


text:
using QMTK
using Compat.Test

@testset "Utils" begin

using QMTK.Consts.Pauli

@test kronprod(sigmax, sigmax, sigmai) == kron(kron(sigmax, sigmax), sigmai)
@test sigmax ⊗ sigmay ⊗ sigmaz == kron(kron(sigmax, sigmay), sigmaz)

h = @kron sigmax[1] ⊗ sigmaz[3] + sigmax[2] ⊗ sigmay[4]
ans = kronprod(σ₁, σ₀, σ₃, σ₀) + kronprod(σ₀, σ₁, σ₀, σ₂)
@test h == ans

h = @kron sigmax[1] * (sigmaz[2] + sigmay[4])
ans = kronprod(σ₁, σ₃, σ₀, σ₀) + kronprod(σ₁, σ₀, σ₀, σ₂)
@test h == ans

h = @kron sigmax[5] * (sigmax[1] * (sigmaz[2] + sigmay[4]))
ans = kronprod(σ₁, σ₃, σ₀, σ₀, σ₁) + kronprod(σ₁, σ₀, σ₀, σ₂, σ₁)
@test h == ans

end


text: twosort(array) = join(minimum(array),"***")

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
# Assuming all SNPs are biallelic.
function make_binary(haplotypes::AbstractVector{T}) where T <: AbstractVector
    n = length(haplotypes)
    n ≡ 0 && return zeros(Bool, 0, 0)
    h1 = haplotypes[1]
    m = length(h1)
    snps = zeros(bool, n, m)
    m ≡ 0 && return snps
    for i in 2:n
        h = haplotypes[i]
        for j in 1:m
            snps[i, j] = h[j] ≢ h1[j]
        end
    end
    snps
end

# Calculate frequencies and make major alleles `true`s.
function prepare_select(haplotypes::Matrix{Bool})
    n, m = size(haplotypes)
    haplotypes = copy(haplotypes)
    fs = vec(sum(haplotypes; dims = 1) / n)
    for j in 1:m
        if fs[j] < 0.5
            haplotypes[:, j] .= .!haplotypes[:, j]
            fs[j] = 1 - fs[j]
        end
    end
    fs, haplotypes
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

@testset "graph.multigraph" begin
    p5 = pathgraph(5)
    @test !hasselfloop(p5)
    @test !hasmultiedge(p5)
    @test issimplegraph(p5)
    looped = immutableplaingraph(4, [(1, 2), (2, 3), (3, 1), (3, 3), (3, 4)])
    @test hasselfloop(looped)
    @test !hasmultiedge(looped)
    @test !issimplegraph(looped)
    multi = immutableplaingraph(4, [(1, 2), (2, 3), (3, 1), (3, 4), (3, 4)])
    @test !hasselfloop(multi)
    @test hasmultiedge(multi)
    @test !issimplegraph(multi)
end # graph.multigraph


text:
# ---
# title: 707. Design Linked List
# id: problem707
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List, Design
# link: <https://leetcode.com/problems/design-linked-list/description/>
# hidden: true
# ---
#
# Design your implementation of the linked list. You can choose to use a singly
# or doubly linked list.
# A node in a singly linked list should have two attributes: `val` and `next`.
# `val` is the value of the current node, and `next` is a pointer/reference to
# the next node.
# If you want to use the doubly linked list, you will need one more attribute
# `prev` to indicate the previous node in the linked list. Assume all nodes in
# the linked list are **0-indexed**.
#
# Implement the `MyLinkedList` class:
#
#   * `MyLinkedList()` Initializes the `MyLinkedList` object.
#   * `int get(int index)` Get the value of the `indexth` node in the linked list. If the index is invalid, return `-1`.
#   * `void addAtHead(int val)` Add a node of value `val` before the first element of the linked list. After the insertion, the new node will be the first node of the linked list.
#   * `void addAtTail(int val)` Append a node of value `val` as the last element of the linked list.
#   * `void addAtIndex(int index, int val)` Add a node of value `val` before the `indexth` node in the linked list. If `index` equals the length of the linked list, the node will be appended to the end of the linked list. If `index` is greater than the length, the node **will not be inserted**.
#   * `void deleteAtIndex(int index)` Delete the `indexth` node in the linked list, if the index is valid.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["MyLinkedList", "addAtHead", "addAtTail", "addAtIndex", "get", "deleteAtIndex", "get"]
#     [[], [1], [3], [1, 2], [1], [1], [1]]
#     **Output**
#     [null, null, null, null, 2, null, 3]
#
#     **Explanation**
#     MyLinkedList myLinkedList = new MyLinkedList();
#     myLinkedList.addAtHead(1);
#     myLinkedList.addAtTail(3);
#     myLinkedList.addAtIndex(1, 2);    // linked list becomes 1->2->3
#     myLinkedList.get(1);              // return 2
#     myLinkedList.deleteAtIndex(1);    // now the linked list is 1->3
#     myLinkedList.get(1);              // return 3
#
#
#
#
# **Constraints:**
#
#   * `0 <= index, val <= 1000`
#   * Please do not use the built-in LinkedList library.
#   * At most `2000` calls will be made to `get`, `addAtHead`, `addAtTail`,  `addAtIndex` and `deleteAtIndex`.
#
#
## @lc code=start
using LeetCode

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


text:
using .Plots
using DSP
using DSP.Windows

"""
    plot(data::SampledSignal; kwargs...)

Plots timeseries from a sample buffer.

# Optional keyword arguments
- `t0=0.0`: start time
- `downsample=:auto`: downsampling factor (integer)
- `pooling=:auto`: pooling mode (`:min`, `:max`, `:mean`, `:minmax`, `nothing` or function)
- other `kwargs` are passed on to `plot`
"""
@recipe function plot(s::SampledSignal; t0=0.0, downsample=:auto, pooling=:auto)
  ticks --> :native
  legend --> ndims(s) > 1 && size(s, 2) > 1
  s1 = samples(s)
  t = domain(s) .+ t0
  if maximum(t) <= 1.0
    t *= 1000.0
    xguide --> "Time (ms)"
  else
    xguide --> "Time (s)"
  end
  n = size(s1, 1)
  if downsample == :auto
    downsample = n >= 10000 ? ceil(Int, n/5000) : 1
    #downsample > 1 && @warn "Downsampling by $downsample"
  end
  if downsample !== nothing && downsample != 1
    if pooling === nothing || (pooling == :auto && isanalytic(s1))
      phase = floor(Int, downsample/2)
      s1 = s1[1+phase:downsample:end,:]
      t = t[1+phase:downsample:end]
    else
      if pooling == :minmax || pooling == :auto
        pooling = orderedextrema
      elseif pooling == :min
        pooling = minimum
      elseif pooling == :max
        pooling = maximum
      elseif pooling == :mean
        pooling = mean
      end
      if ndims(s) == 1
        s1 = collect(Iterators.flatten(pooling.(Periodograms.arraysplit(s1, downsample, 0))))
        m = fld(n, downsample)
        rep, rem = fldmod(size(s1, 1), m)
        rem == 0 || throw(ArgumentError("Pooling function must sample the signal evenly"))
        t = t[1:downsample:end]
        rep > 1 && (t = repeat(t, inner=rep))
      else
        y = collect(Iterators.flatten(pooling.(Periodograms.arraysplit(s1[:,1], downsample, 0))))
        if size(s1, 2) > 1
          y1 = zeros(length(y), size(s1, 2))
          y1[:,1] .= y
          y = y1
        end
        for j in 2:size(s1, 2)
          y[:,j] .= collect(Iterators.flatten(pooling.(Periodograms.arraysplit(s1[:,j], downsample, 0))))
        end
        s1 = y
        m = fld(n, downsample)
        rep, rem = fldmod(size(s1, 1), m)
        rem == 0 || throw(ArgumentError("Pooling function must sample the signal evenly"))
        t = t[1:downsample:end]
        rep > 1 && (t = repeat(t, inner=rep))
      end
    end
  end
  @series begin
    seriestype := :line
    t[1:size(s1,1)], s1
  end
end

"""
    plot(data::TFR; kwargs...)

Displays time frequency representation.

# Optional keyword arguments
- `t0=0.0`: start time
- `cscale=:log`: color scale type (`:log`, `:linear`)
- `crange=:auto`: color scale dynamic range (`:auto`, number)
- other `kwargs` are passed on to `plot`
"""
@recipe function plot(tfr::Periodograms.TFR; t0=0.0, cscale=:log, crange=:auto)
  ticks --> :native
  t = time(tfr) .+ t0
  f = freq(tfr)
  z = power(tfr)
  if maximum(t) <= 1.0
    t *= 1000.0
    xguide --> "Time (ms)"
  else
    xguide --> "Time (s)"
  end
  if maximum(f) >= 5000.0
    f /= 1000.0
    yguide --> "Frequency (kHz)"
  else
    yguide --> "Frequency (Hz)"
  end
  z = abs.(z)
  if cscale == :log
    z = pow2db.(z)
    cmax = ceil(maximum(z)/5)*5
    crange == :auto && (crange = 50)
    cmin = max(cmax-crange, floor(minimum(z)/5)*5)
    clims --> (cmin, cmax)
  elseif crange != :auto
    cmax = maximum(z)
    cmin = max(cmax-crange, minimum(z))
    clims --> (cmin, cmax)
  end
  @series begin
    seriestype := :heatmap
    t, f, z
  end
end

"""
    psd(data; kwargs...)

Plots the power spectral density of data.

# Optional keyword arguments
- `fs=1.0`: derived from the data if a `SampledSignal` is provided as input
- `nfft=512`: size of FFT window
- `noverlap=nfft÷2`: window overlap size
- `window=hamming(nfft)`: accepts any window from [`DSP.jl`](https://juliadsp.github.io/DSP.jl/stable/windows/)
- `xscale=:auto`: one of `:auto`, `:identity` or `:log10`
- `yrange=50`: y-scale dB range for automatic scaling
- other `kwargs` are passed on to `plot`
"""
@userplot PSD
@recipe function plot(s::PSD; fs=1.0, nfft=512, noverlap=div(nfft,2),
                      window=hamming(nfft), xscale=:auto, yrange=50)
  if length(s.args) != 1 || !(s.args[1] isa AbstractArray)
    throw(ArgumentError("psd should be provided timeseries data"))
  end
  s.args[1] isa SampledSignal && (fs = framerate(s.args[1]))
  s = samples(s.args[1])
  nfft = nextfastfft(nfft)
  while nfft > size(s, 1)
    nfft = nfft ÷ 2
  end
  ticks --> :native
  legend --> ndims(s) > 1 && size(s, 2) > 1
  yguide --> "Power spectral density (dB/Hz)"
  xscale == :auto && (xscale := :identity)
  for j = 1:size(s, 2)
    p = welch_pgram(s[:,j], nfft, noverlap; fs=inHz(fs), window=window)
    f = freq(p)
    if maximum(f) > 10000 && xscale == :auto
      xguide --> "Frequency (kHz)"
      f /= 1000.0
    else
      xguide --> "Frequency (Hz)"
    end
    if xscale == :log10
      xlims --> (f[2], f[end])
    end
    y = pow2db.(power(p))
    ymax = ceil(Int, maximum(y)/5)*5
    ymin = max(ymax-yrange, floor(Int, minimum(y)/5)*5)
    ylims --> (ymin, ymax)
    @series begin
      seriestype := :line
      f, y
    end
  end
end

"""
    specgram(data; kwargs...)

Plots a spectrogram of the data.

# Optional keyword arguments
- `fs=1.0`: derived from the data if a `SampledSignal` is provided as input
- `nfft=256`: size of FFT window
- `noverlap=nfft÷2`: window overlap size
- `window=hamming(nfft)`: accepts any window from [`DSP.jl`](https://juliadsp.github.io/DSP.jl/stable/windows/)
- `t0=0.0`: start time
- `downsample=:auto`: downsampling factor (integer) for time axis
- `pooling=:mean`: pooling mode (`:min`, `:max`, `:mean`, `nothing` or function)
- `crange=50`: color scale dB range for automatic scaling
- other `kwargs` are passed on to `plot`
"""
@userplot Specgram
@recipe function plot(s::Specgram; fs=1.0, nfft=256, noverlap=div(nfft,2), crange=50,
                      window=hamming(nfft), t0=0.0, downsample=:auto, pooling=:mean)
  if length(s.args) != 1 || !(s.args[1] isa AbstractArray)
    throw(ArgumentError("specgram should be provided timeseries data"))
  end
  s.args[1] isa SampledSignal && (fs = framerate(s.args[1]))
  s1 = samples(s.args[1])
  if ndims(s1) > 1 && size(s1, 2) > 1
    throw(ArgumentError("specgram does not support multichannel data"))
  end
  nfft = nextfastfft(nfft)
  while nfft > size(s1, 1)
    nfft = nfft ÷ 2
  end
  ticks --> :native
  p = spectrogram(s1[:,1], nfft, noverlap; fs=inHz(fs), window=window)
  t = t0 .+ time(p)
  f = freq(p)
  z = power(p)
  if downsample == :auto
    downsample = length(t) > 10000 ? ceil(Int, length(t)/5000) : 1
    #downsample > 1 && @warn "Downsampling by $downsample"
  end
  if downsample !== nothing && downsample != 1
    if pooling === nothing
      phase = floor(Int, downsample/2)
      t = t[1+phase:downsample:end]
      z = z[:,1+phase:downsample:end]
    else
      if pooling == :min
        pooling = minimum
      elseif pooling == :max
        pooling = maximum
      elseif pooling == :mean
        pooling = mean
      end
      y1 = collect(Iterators.flatten(pooling.(Periodograms.arraysplit(z[1,:], downsample, 0))))
      y = zeros(size(z,1), length(y1))
      y[1,:] .= y1
      for j in 2:size(z,1)
        y[j,:] .= collect(Iterators.flatten(pooling.(Periodograms.arraysplit(z[j,:], downsample, 0))))
      end
      z = y
      t = t[1:downsample:end][1:size(z,2)]
    end
  end
  if isanalytic(s1)
    ndx = sortperm(f)
    f = f[ndx]
    z = z[ndx,:]
  end
  tunit = "s"
  funit = "Hz"
  if maximum(f) >= 10000
    f /= 1000.0
    funit = "kHz"
  end
  if maximum(t) <= 1.0
    t *= 1000.0
    tunit = "ms"
  end
  xguide --> "Time ("*tunit*")"
  yguide --> "Frequency ("*funit*")"
  z = pow2db.(z)
  cmax = ceil(maximum(z)/5)*5
  cmin = max(cmax-crange, floor(minimum(z)/5)*5)
  clims --> (cmin, cmax)
  @series begin
    seriestype := :heatmap
    t, f, z
  end
end

"""
    plotfreqresp(filter; kwargs...)
    plotfreqresp(num; kwargs...)
    plotfreqresp(num, den; kwargs...)

Plots frequency response of a digital filter.

# Optional keyword arguments
- `fs=1.0`: sampling frequency
- `nfreq=256`: number of frequency points to evaluate filter response at
- `xscale=:auto`: one of `:auto`, `:identity` or `:log10`
- other `kwargs` are passed on to `plot`
"""
@userplot PlotFreqResp
@recipe function plot(s::PlotFreqResp; fs=1.0, nfreq=256, xscale=:auto)
  if length(s.args) ∉ (1, 2)
    throw(ArgumentError("plotfreqresp should be provided filter details"))
  end
  if s.args[1] isa FilterCoefficients
    filt = s.args[1]
  else
    num = s.args[1] isa AbstractArray ? s.args[1] : [s.args[1]]
    if length(s.args) == 2
      den = s.args[2] isa AbstractArray ? s.args[2] : [s.args[2]]
      filt = PolynomialRatio(num, den)
    else
      filt = PolynomialRatio(num, [1.0])
    end
  end
  fs = inHz(fs)
  f = collect(range(0, stop=fs/2, length=nfreq))
  z = Filters.freqresp(filt, 2π * f/fs)
  funit = "Hz"
  if xscale == :auto
    xscale := :identity
    if maximum(f) >= 10000
      f /= 1000.0
      funit = "kHz"
    end
  end
  if xscale == :log10
    f = f[2:end]
    z = z[2:end]
  end
  layout := @layout [mag; phase]
  @series begin
    seriestype := :line
    subplot := 1
    xguide --> "Frequency ("*funit*")"
    yguide --> "Response (dB)"
    legend --> false
    f, amp2db.(abs.(z))
  end
  @series begin
    seriestype := :line
    subplot := 2
    xguide --> "Frequency ("*funit*")"
    yguide --> "Phase (radians)"
    legend --> false
    f, angle.(z)
  end
end

# helper functions

orderedextrema(x) = argmin(x) < argmax(x) ? (minimum(x), maximum(x)) : (maximum(x), minimum(x))


text:
type JointOrderStatistic{D,S<:MultivariateDistribution} <: AbstractJointOrderStatistic{D}
    sequence::S
    orders::Vector{Int64}
end

function JointOrderStatistic{D,S}(sequence::IIDRandomSequence{D,S},orders::Vector{Int64})
    return JointOrderStatistic{D,IIDRandomSequence{D,S}}(sequence,orders)
end

parentsequence(X::JointOrderStatistic) = X.sequence
orders(X::JointOrderStatistic) = X.orders

length(X::JointOrderStatistic) = length(orders(X))

function _rand!{T<:Real}(X::JointOrderStatistic, x::AbstractVector{T})
  y = sort(rand(parentsequence(X)))
  x = copy(y[orders(X)])
  return x
end

typealias IIDJointOrderStatistic{D<:ValueSupport,T<:UnivariateDistribution} JointOrderStatistic{D,IIDRandomSequence{D,T}}
typealias IIDContinuousJointOrderStatistic{T<:ContinuousUnivariateDistribution} IIDJointOrderStatistic{Continuous,T}
typealias IIDDiscreteJointOrderStatistic{T<:DiscreteUnivariateDistribution} IIDJointOrderStatistic{Discrete,T}

function _logpdf{T<:Real,S<:ContinuousUnivariateDistribution}(X::IIDContinuousJointOrderStatistic{S}, x::AbstractVector{T})
  N,ks = length(parentsequence(X)),orders(X)
  logps,Ps = logpdf(X.sequence.d,x),cdf(X.sequence.d,x)
  ΔPs,Δqs = diff([0.;Ps;1.]),diff([0;ks;ks[end]+1])-1
  return lfact(N)+sum(logps)+sum(Δqs.*log(ΔPs))- sum(lfact(Δqs))
end

function mean(X::IIDJointOrderStatistic)
  return map(mean,[ScalarOrderStatistic(parentsequence(X),order) for order in orders(X)])
end


text:
"""
    datasetpath(name)

Return the folder that dataset `name` is stored.

If it hasn't been downloaded yet, you will be asked if you want to
download it. See [`Datasets.DATASETS`](#) for a list of available datasets.
"""
function datasetpath(name)
    datadeppath = @datadep_str "fastai-$name"
    return Path(joinpath(datadeppath, name))
end

"""
    loadtaskdata(dir, Task)
    loadtaskdata(dir, method::LearningMethod{Task})

Load a task data container for `LearningTask` `Task` stored in `dir`
in a canonical format.
"""
function loadtaskdata(dir, ::DLPipelines.LearningMethod{T}) where T
    return loadtaskdata(dir, T)
end

"""
    loadtaskdata(dir, ImageClassificationTask; [split = false])

Load a data container for `ImageClassificationTask` with observations
`(input = image, target = class)`.

If `split` is `true`, returns a tuple of the data containers split by
the name of the grandparent folder.

`dir` should contain the data in the following canonical format:

- dir
    - split (e.g. "train", "valid"...)
        - class (e.g. "cat", "dog"...)
            - image32434.{jpg/png/...}
            - ...
        - ...
    - ...
"""
function loadtaskdata(
        dir,
        Task::Type{FastAI.ImageClassificationTask};
        split=false,
        filterparent= (!=("test")),
        kwargs...)
    data = filterobs(FileDataset(dir)) do path
        isimagefile(path) && filterparent(filename(parent(path)))
    end
    if split
        datas = groupobs(data) do path
            filename(parent(parent(path)))
        end
        return Tuple(mapobs(
            (input = loadfile, target = path -> filename(parent(path))),
            data
        ) for data in datas)
    else
        return mapobs(
            (input = loadfile, target = path -> filename(parent(path))),
            data
        )
    end
end

"""
    getclasses(name)

Get the list of classes for classification dataset `name`.
"""
function getclassesclassification(dir)
    data = mapobs(filterobs(isimagefile, FileDataset(dir))) do path
        return filename(parent(path))
    end
    return unique(collect(eachobsparallel(data, useprimary=true, buffered=false)))
end
getclassesclassification(name::String) = getclassesclassification(datasetpath(name))

"""
    getclasses(name)

Get the list of classes for classification dataset `name`.
"""
function getclassessegmentation(dir::AbstractPath)
    classes = readlines(open(joinpath(dir, "codes.txt")))
    return classes
end
getclassessegmentation(name::String) = getclassessegmentation(datasetpath(name))

"""
    loadtaskdata(dir, ImageSegmentationTask; [split = false])

Load a data container for `ImageSegmentationTask` with observations
`(input = image, target = mask)`.

If `split` is `true`, returns a tuple of the data containers split by
the name of the grandparent folder.

"""
function loadtaskdata(
        dir,
        Task::Type{FastAI.ImageSegmentationTask};
        split=false,
        kwargs...)
    imagedata = mapobs(loadfile, filterobs(isimagefile, FileDataset(joinpath(dir, "images"))))
    maskdata = mapobs(maskfromimage ∘ loadfile, filterobs(isimagefile, FileDataset(joinpath(dir, "labels"))))
    return mapobs((input = obs -> obs[1], target = obs -> obs[2]), (imagedata, maskdata))
end

maskfromimage(a::AbstractArray{<:Gray{T}}) where T = maskfromimage(reinterpret(T, a))
maskfromimage(a::AbstractArray{<:Normed{T}}) where T = maskfromimage(reinterpret(T, a))
function maskfromimage(a::AbstractArray{I}) where {I<:Integer}
    return a .+ one(I)
end


text:
TypeObject = Union{DataType, Union, UnionAll}

"""representing the in-matching object in pattern compile time
"""
struct Target{IsComplex}
    repr::Any
    type::Ref{TypeObject}
end

@nospecialize

function target_method end

function Base.getproperty(target::Target, sym::Symbol)
    target_method(target, Val(sym))
end

function target_method(target::Target, ::Val{:type_narrow!})
    function (ty::TypeObject)
        tyref = getfield(target, :type)
        if tyref[] <: ty
        else
            tyref[] = ty
        end
    end
end

function target_method(target::Target, ::Val{:type})
    getfield(target, :type)[]
end

function target_method(target::Target, ::Val{:repr})
    getfield(target, :repr)
end

function target_method(target::Target{IsC}, ::Val{:with_repr}) where IsC
    function ap(repr::Any)
        Target{IsC}(repr,  getfield(target, :type))
    end
    function ap(repr::Any, ::Val{IsC′}) where IsC′
        Target{IsC′}(repr, getfield(target, :type))
    end
    ap
end

function target_method(target::Target{IsC}, ::Val{:with_type}) where IsC
    function ap(ty::TypeObject)
        Target{IsC}(target.repr, Ref{TypeObject}(ty))
    end
    function ap(ty::TypeObject, ::Val{IsC′}) where IsC′
        Target{IsC′}(target.repr, Ref{TypeObject}(ty))
    end
    ap
end

function target_method(target::Target{IsC}, ::Val{:clone}) where IsC
    Target{IsC}(target.repr, Ref{TypeObject}(target.type))
end

@specialize


text:
############################################################################################
"""
$(TYPEDEF)
Abstract super type for Supported Automatic Differentiation backends.
"""
abstract type AutomaticDifferentiationMethod end
struct ADForward <: AutomaticDifferentiationMethod end
struct ADZygote <: AutomaticDifferentiationMethod end
struct ADReverse <: AutomaticDifferentiationMethod end
struct ADReverseUntaped <: AutomaticDifferentiationMethod end

############################################################################################
#Helper function for AD packages
function checkfinite(θ::ReverseDiff.TrackedArray{T}, max_val::R=max_val) where {T,R<:Real}
    @inbounds @simd for iter in eachindex(θ)
        if !checkfinite(θ[iter], max_val=max_val)
            return false
        end
    end
    return true
end

"""
$(SIGNATURES)
Initiate DiffResults.MutableDiffResult struct. Not exported.

# Examples
```julia
```

"""
function _diffresults_buffer(θᵤ::AbstractVector{T}) where {T<:Real}
    #NOTE: Adjusted from: https://github.com/tpapp/LogDensityProblems.jl/blob/master/src/DiffResults_helpers.jl
    S = T <: Real ? T : Float64
    return DiffResults.MutableDiffResult(zero(S), (Vector{S}(undef, size(θᵤ, 1)),))
end
############################################################################################
struct AutomaticDiffTune{M<:AutomaticDifferentiationMethod,C} <: AbstractDifferentiableTune
    "Automatic Differentiation (AD) backend."
    backend::M
    "Chunck size configuration for AD backend."
    config::C
    function AutomaticDiffTune(backend::M, config::C
    ) where {M<:AutomaticDifferentiationMethod,C}
        return new{M,C}(backend, config)
    end
end
function AutomaticDiffTune(backend::M, objective::Objective) where {M<:AutomaticDifferentiationMethod}
    return AutomaticDiffTune(
    backend, _config(backend, objective, unconstrain_flatten(objective.model, objective.tagged))
    )
end

############################################################################################
function update(
    tune::AutomaticDiffTune,
    objective::Objective)
    return AutomaticDiffTune(tune.backend, objective)
end

############################################################################################
function _log_density(
    objective::Objective, tune::AutomaticDiffTune, θᵤ::AbstractVector{T}
) where {T<:Real}
    return objective(θᵤ)
end

############################################################################################
#Make Symbol intialization easier
function AutomaticDiffTune(backend::Symbol, objective::Objective)
    return AutomaticDiffTune(Val(backend), objective)
end

############################################################################################
#!NOTE: All those functions need to be dispatched if new AD framework is included
"""
$(SIGNATURES)
Write config file for AD wrapper.

# Examples
```julia
```

"""
function _config(
    differentiation::ADForward, objective::Objective, θᵤ::AbstractVector{R}
) where {R<:Real}
    return ForwardDiff.GradientConfig(objective, θᵤ)
end
function _config(
    differentiation::ADZygote, objective::Objective, θᵤ::AbstractVector{R}
) where {R<:Real}
    return nothing
end
function _config(
    differentiation::ADReverse, objective::Objective, θᵤ::AbstractVector{R}
) where {R<:Real}
    return ReverseDiff.compile(ReverseDiff.GradientTape(objective, θᵤ))
end
function _config(
    differentiation::ADReverseUntaped, objective::Objective, θᵤ::AbstractVector{R}
) where {R<:Real}
    return nothing
end

############################################################################################
function AutomaticDiffTune(
    backend::Val{:ForwardDiff},
    objective::Objective,
    config::C=_config(ADForward(), objective, unconstrain_flatten(objective.model, objective.tagged)),
) where {C}
    return AutomaticDiffTune(ADForward(), config)
end

function AutomaticDiffTune(
    backend::Val{:Zygote},
    objective::Objective,
    config::C=_config(ADZygote(), objective, unconstrain_flatten(objective.model, objective.tagged)),
) where {C}
    return AutomaticDiffTune(ADZygote(), config)
end

function AutomaticDiffTune(
    backend::Val{:ReverseDiff},
    objective::Objective,
    config::C=_config(ADReverse(), objective, unconstrain_flatten(objective.model, objective.tagged)),
) where {C}
    return AutomaticDiffTune(ADReverse(), config)
end

function AutomaticDiffTune(
    backend::Val{:ReverseDiffUntaped},
    objective::Objective,
    config::C=_config(ADReverseUntaped(), objective, unconstrain_flatten(objective.model, objective.tagged)),
) where {C}
    return AutomaticDiffTune(ADReverseUntaped(), config)
end

############################################################################################
function _log_density_and_gradient(
    objective::Objective, tune::AutomaticDiffTune{ADForward}, θᵤ::AbstractVector{T}
) where {T<:Real}
    buffer = _diffresults_buffer(θᵤ)
    result = ForwardDiff.gradient!(buffer, objective, θᵤ, tune.config)
    return DiffResults.value(result), DiffResults.gradient(result)
end
function _log_density_and_gradient(
    objective::Objective, tune::AutomaticDiffTune{ADZygote}, θᵤ::AbstractVector{T}
) where {T<:Real}
    _val, back = Zygote.pullback(objective, θᵤ)
    return T(_val), first(back(Zygote.sensitivity(_val)))
end
function _log_density_and_gradient(
    objective::Objective, tune::AutomaticDiffTune{ADReverse}, θᵤ::AbstractVector{T}
) where {T<:Real}
    buffer = _diffresults_buffer(θᵤ)
    result = ReverseDiff.gradient!(buffer, tune.config, θᵤ)
    return DiffResults.value(result), DiffResults.gradient(result)
end
function _log_density_and_gradient(
    objective::Objective, tune::AutomaticDiffTune{ADReverseUntaped}, θᵤ::AbstractVector{T}
) where {T<:Real}
    buffer = _diffresults_buffer(θᵤ)
    result = ReverseDiff.gradient!(buffer, objective, θᵤ)
    return DiffResults.value(result), DiffResults.gradient(result)
end

############################################################################################
# Export
export AutomaticDiffTune


text:
using Test, MPI
include(joinpath("..","..","testhelpers.jl"))

include("iterativesolvers.jl")

@testset "Linear Solvers Poisson" begin
    tests = [
        (1, "columnwiselu.jl"),
        (1, "poisson.jl"),
        (1, "bandedsystem.jl"),
        (1, "cg.jl"),
    ]
    runmpi(tests, @__FILE__)
end


text:
"""
`RandomODEProblem` - initialize the problem data structure

Defaults to the zero potential; in this case the ODE becomes

-Δu = f

on (a, b) with the specified boundary conditions.
"""
function RandomODEProblem(a, b, nₓ, m, σ, s, left_bc, right_bc, f; potential = x -> 0.)
    Δx = (b - a)/(nₓ + 1)
    x = [a + i*Δx for i ∈ 1:nₓ]
    mean = m.(x)
    V = potential.(x)
    A = spzeros(nₓ, nₓ)
    rhs = f.(x)
    u = zeros(nₓ)
    return RandomODEProblem(a, b, nₓ, Δx, x, mean, σ, s, V, left_bc, right_bc, A, f, rhs, u)
end


text:
module LoweredCodeUtils

using Core: SimpleVector, CodeInfo
using Base.Meta: isexpr
using JuliaInterpreter
using JuliaInterpreter: SSAValue, SlotNumber, Frame
using JuliaInterpreter: @lookup, moduleof, pc_expr, step_expr!, is_global_ref, whichtt,
                        next_until!, finish_and_return!, nstatements, codelocation

export signature, methoddef!, methoddefs!, bodymethod

"""
    iscallto(stmt, name)

Returns `true` is `stmt` is a call expression to `name`.
"""
function iscallto(stmt, name)
    if isexpr(stmt, :call)
        a = stmt.args[1]
        a == name && return true
        return is_global_ref(a, Core, :_apply) && stmt.args[2] == name
    end
    return false
end

"""
    getcallee(stmt)

Returns the function (or Symbol) being called in a :call expression.
"""
function getcallee(stmt)
    if isexpr(stmt, :call)
        a = stmt.args[1]
        is_global_ref(a, Core, :_apply) && return stmt.args[2]
        return a
    end
    error(stmt, " is not a call expression")
end

ismethod(frame::Frame) = ismethod(pc_expr(frame))
ismethod(stmt)  = isexpr(stmt, :method)
ismethod1(stmt) = isexpr(stmt, :method, 1)
ismethod3(stmt) = isexpr(stmt, :method, 3)

function methodname(name)
    isa(name, Symbol) && return name
    if isa(name, CodeInfo) && length(name.code) == 1 && isexpr(name.code[1], :return) && ismethod1(name.code[1].args[1])
        return name.code[1].args[1].args[1]
    end
    error("unhandled name type ", typeof(name))
end

"""
    nextpc = next_or_nothing(frame, pc)
    nextpc = next_or_nothing!(frame)

Advance the program counter without executing the corresponding line.
If `frame` is finished, `nextpc` will be `nothing`.
"""
next_or_nothing(frame, pc) = pc < nstatements(frame.framecode) ? pc+1 : nothing
function next_or_nothing!(frame)
    pc = frame.pc
    if pc < nstatements(frame.framecode)
        frame.pc = pc = pc + 1
        return pc
    end
    return nothing
end

"""
    nextpc = skip_until(predicate, frame, pc)
    nextpc = skip_until!(predicate, frame)

Advance the program counter until `predicate(stmt)` return `true`.
"""
function skip_until(predicate, frame, pc)
    stmt = pc_expr(frame, pc)
    while !predicate(stmt)
        pc = next_or_nothing(frame, pc)
        pc === nothing && return nothing
        stmt = pc_expr(frame, pc)
    end
    return pc
end
function skip_until!(predicate, frame)
    pc = frame.pc
    stmt = pc_expr(frame, pc)
    while !predicate(stmt)
        pc = next_or_nothing!(frame)
        pc === nothing && return nothing
        stmt = pc_expr(frame, pc)
    end
    return pc
end

"""
    sig = signature(sigsv::SimpleVector)

Compute a method signature from a suitable `SimpleVector`: `sigsv[1]` holds the signature
and `sigsv[2]` the `TypeVar`s.

# Example:

For `f(x::AbstractArray{T}) where T`, the corresponding `sigsv` is constructed as

    T = TypeVar(:T)
    sig1 = Core.svec(typeof(f), AbstractArray{T})
    sig2 = Core.svec(T)
    sigsv = Core.svec(sig1, sig2)
    sig = signature(sigsv)
"""
function signature(sigsv::SimpleVector)
    sigp, sigtv = sigsv
    sig = Tuple{sigp...}
    for i = length(sigtv):-1:1
        sig = UnionAll(sigtv[i], sig)
    end
    return sig
end

"""
    sigt, lastpc = signature(recurse, frame, pc)
    sigt, lastpc = signature(frame, pc)

Compute the signature-type `sigt` of a method whose definition in `frame` starts at `pc`.
Generally, `pc` should point to the `Expr(:method, methname)` statement, in which case
`lastpc` is the final statement number in `frame` that is part of the signature
(i.e, the line above the 3-argument `:method` expression).
Alternatively, `pc` can point to the 3-argument `:method` expression,
as long as all the relevant SSAValues have been assigned.
In this case, `lastpc == pc`.

If no 3-argument `:method` expression is found, `sigt` will be `nothing`.
"""
function signature(@nospecialize(recurse), frame::Frame, @nospecialize(stmt), pc)
    mod = moduleof(frame)
    lastpc = frame.pc = pc
    while !isexpr(stmt, :method, 3)  # wait for the 3-arg version
        if isexpr(stmt, :thunk) && isanonymous_typedef(stmt.args[1])
            lastpc = pc = define_anonymous(recurse, frame, stmt)
        elseif isexpr(stmt, :call) && JuliaInterpreter.is_quotenode(stmt.args[1], Core.Typeof) &&
               (sym = stmt.args[2]; isa(sym, Symbol) && !isdefined(mod, sym))
            return nothing, pc
        else
            lastpc = pc
            pc = step_expr!(recurse, frame, stmt, true)
            pc === nothing && return nothing, lastpc
        end
        stmt = pc_expr(frame, pc)
    end
    sigsv = @lookup(frame, stmt.args[2])::SimpleVector
    sigt = signature(sigsv)
    return sigt, lastpc
end
signature(@nospecialize(recurse), frame::Frame, pc) = signature(recurse, frame, pc_expr(frame, pc), pc)
signature(frame::Frame, pc) = signature(finish_and_return!, frame, pc)

function minid(node, stmts, id)
    if isa(node, SSAValue)
        id = min(id, node.id)
        stmt = stmts[node.id]
        return minid(stmt, stmts, id)
    elseif isa(node, Expr)
        for a in node.args
            id = minid(a, stmts, id)
        end
    end
    return id
end

function signature_top(frame, stmt, pc)
    @assert ismethod3(stmt)
    return minid(stmt.args[2], frame.framecode.src.code, pc)
end

##
## Detecting anonymous functions. These start with a :thunk expr and have a characteristic CodeInfo
##
function isanonymous_typedef(code::CodeInfo)
    length(code.code) >= 4 || return false
    stmt = code.code[end-1]
    isexpr(stmt, :struct_type) || return false
    name = stmt.args[1]::Symbol
    return startswith(String(name), "#")
end

function define_anonymous(@nospecialize(recurse), frame, stmt)
    while !isexpr(stmt, :method)
        pc = step_expr!(recurse, frame, stmt, true)
        stmt = pc_expr(frame, pc)
    end
    return step_expr!(recurse, frame, stmt, true)  # also define the method
end

##
## Deal with gensymmed names, https://github.com/JuliaLang/julia/issues/30908
## This handles the situation in which a method is created using a different gensym
## than when this method was last lowered & evaled.
##

function get_parentname(name)
    isa(name, Expr) && return name
    name = methodname(name)
    isa(name, Symbol) || error("unhandled name type ", typeof(name))
    namestring = String(name)
    if namestring[1] == '#'
        pidx = 1
        idx = nextind(namestring, pidx)
        while namestring[idx] != '#'
            pidx = idx
            idx = nextind(namestring, idx)
        end
        parentname = Symbol(namestring[2:pidx])
    else
        parentname = name
    end
    return parentname
end

"""
    pctop, isgen = find_corrected_name(frame, pc, name, parentname)

Scans forward from `pc` in `frame` until a method is found that calls `name`.
`pctop` points to the beginning of that method's signature.
`isgen` is true if `name` corresponds to sa GeneratedFunctionStub.

Alternatively, this returns `nothing` if `pc` does not appear to point to either
a keyword or generated method.
"""
function find_corrected_name(frame, pc, name, parentname)
    stmt = pc_expr(frame, pc)
    while true
        while !ismethod3(stmt)
            pc = next_or_nothing(frame, pc)
            pc === nothing && return nothing
            stmt = pc_expr(frame, pc)
        end
        body = stmt.args[3]
        if stmt.args[1] != name && stmt.args[1] != parentname
            # This might be the GeneratedFunctionStub for a @generated method
            bodystmt = body.code[1]
            if isexpr(bodystmt, :meta) && bodystmt.args[1] == :generated
                return signature_top(frame, stmt, pc), true
            end
        elseif length(body.code) > 1
            bodystmt = body.code[end-1]  # the line before the final return
            iscallto(bodystmt, name) && return signature_top(frame, stmt, pc), false
        end
        pc = next_or_nothing(frame, pc)
        pc === nothing && return nothing
        stmt = pc_expr(frame, pc)
    end
end

"""
    replacename!(stmts, oldname=>newname)

Replace a Symbol `oldname` with `newname` in `stmts`.
"""
function replacename!(ex::Expr, pr)
    replacename!(ex.args, pr)
    return ex
end

function replacename!(args::Vector{Any}, pr)
    oldname, newname = pr
    for i = 1:length(args)
        a = args[i]
        if isa(a, Expr)
            replacename!(a, pr)
        elseif isa(a, CodeInfo)
            replacename!(a.code, pr)
        elseif isa(a, QuoteNode) && a.value == oldname
            args[i] = QuoteNode(newname)
        elseif a == oldname
            args[i] = newname
        end
    end
    return args
end

function correct_name!(@nospecialize(recurse), frame, pc, name, parentname)
    # Get the correct name (the one that's actively running)
    nameinfo = find_corrected_name(frame, pc, name, parentname)
    if nameinfo === nothing
        pc = skip_until(stmt->isexpr(stmt, :method, 3), frame, pc)
        pc = next_or_nothing(frame, pc)
        return name, pc
    end
    pctop, isgen = nameinfo
    sigtparent, lastpcparent = signature(recurse, frame, pctop)
    sigtparent === nothing && return name, pc
    methparent = whichtt(sigtparent)
    methparent === nothing && return name, pc  # caller isn't defined, no correction is needed
    if isgen
        cname = nameof(methparent.generator.gen)
    else
        bodyparent = Base.uncompressed_ast(methparent)
        bodystmt = bodyparent.code[end-1]
        @assert isexpr(bodystmt, :call)
        ref = getcallee(bodystmt)
        isa(ref, GlobalRef) || @show ref typeof(ref)
        @assert isa(ref, GlobalRef)
        @assert ref.mod == moduleof(frame)
        cname = ref.name
    end
    # Swap in the correct name
    if name != cname
        replacename!(frame.framecode.src.code, name=>cname)
    end
    stmt = pc_expr(frame, lastpcparent)
    while !ismethod(stmt)
        lastpcparent = next_or_nothing(frame, lastpcparent)
        lastpcparent === nothing && return name, lastpcparent
        stmt = pc_expr(frame, lastpcparent)
    end
    name = stmt.args[1]
    return name, pc
end

"""
    ret = methoddef!(recurse, signatures, frame; define=true)
    ret = methoddef!(signatures, frame; define=true)

Compute the signature of a method definition. `frame.pc` should point to a
`:method` expression. Upon exit, the new signature will be added to `signatures`.

There are several possible return values:

    pc, pc3 = ret

is the typical return. `pc` will point to the next statement to be executed, or be `nothing`
if there are no further statements in `frame`. `pc3` will point to the 3-argument `:method`
expression.

Alternatively,

    pc = ret

occurs for "empty method" expressions, e.g., `:(function foo end)`. `pc` will be `nothing`.

By default the method will be defined (evaluated). You can prevent this by setting `define=false`.
This is recommended if you are simply extracting signatures from code that has already been evaluated.

An important point is that `methoddef!` may, in some circumstances, change the names
of methods in `frame`.  The issues are described in https://github.com/JuliaLang/julia/issues/30908.
`methoddef!` will try to replace names with the ones that are currently active.
"""
function methoddef!(@nospecialize(recurse), signatures, frame::Frame, @nospecialize(stmt), pc::Int; define=true)
    framecode = frame.framecode
    if ismethod3(stmt)
        pc3 = pc
        sigt, pc = signature(recurse, frame, stmt, pc)
        if sigt === nothing && define
            step_expr!(recurse, frame, stmt, true)
        end
        sigt, pc = signature(recurse, frame, stmt, pc)
        meth = whichtt(sigt)
        if meth === nothing && define
            step_expr!(recurse, frame, stmt, true)
            meth = whichtt(sigt)
        end
        if isa(meth, Method)
            push!(signatures, meth.sig)
        elseif stmt.args[1] == false
            # If it's anonymous and not defined, define it
            pc = step_expr!(recurse, frame, stmt, true)
            meth = whichtt(sigt)
            isa(meth, Method) && push!(signatures, meth.sig)
            return pc, pc3
        else
            # guard against busted lookup, e.g., https://github.com/JuliaLang/julia/issues/31112
            code = framecode.src
            codeloc = codelocation(code, pc)
            loc = code.linetable[codeloc]
            ft = Base.unwrap_unionall(Base.unwrap_unionall(sigt).parameters[1])
            if !startswith(String(ft.name.name), "##")
                @warn "file $(loc.file), line $(loc.line): no method found for $sigt"
            end
        end
        frame.pc = pc
        return ( define ? step_expr!(recurse, frame, stmt, true) : next_or_nothing!(frame) ), pc3
    end
    ismethod1(stmt) || error("expected method opening, got ", stmt)
    name = stmt.args[1]
    if isa(name, Bool)
        error("not valid for anonymous methods")
    end
    parentname = get_parentname(name)  # e.g., name = #foo#7 and parentname = foo
    pcinc = 1
    nextstmt = pc_expr(frame, pc+pcinc)
    while ismethod1(nextstmt) || isexpr(nextstmt, :global)
        if ismethod1(nextstmt)
            name = nextstmt.args[1]
        end
        pcinc += 1
        nextstmt = pc_expr(frame, pc+pcinc)
    end
    if !define && String(name)[1] == '#'
        # We will have to correct the name.
        # We can only correct one at a time, so work backwards from a non-gensymmed name
        # (https://github.com/timholy/Revise.jl/issues/290)
        pc0 = pc
        idx1 = findall(ismethod1, frame.framecode.src.code)
        idx1 = idx1[idx1 .>= pc]
        hashhash = map(idx->startswith(String(pc_expr(frame, idx).args[1]), '#'), idx1)
        idx1 = idx1[hashhash]
        i = length(idx1)
        while i > 1
            frame.pc = idx1[i]
            methoddef!(recurse, [], frame, frame.pc; define=define)
            i -= 1
        end
        frame.pc = pc0
    end
    if name != parentname && !define
        name, endpc = correct_name!(recurse, frame, pc, name, parentname)
    end
    while true  # methods containing inner methods may need multiple trips through this loop
        sigt, pc = signature(recurse, frame, stmt, pc)
        stmt = pc_expr(frame, pc)
        while !isexpr(stmt, :method, 3)
            pc = next_or_nothing(frame, pc)  # this should not check define, we've probably already done this once
            pc === nothing && return nothing   # this was just `function foo end`, signal "no def"
            stmt = pc_expr(frame, pc)
        end
        pc3 = pc
        name3 = stmt.args[1]
        sigt === nothing && (error("expected a signature"); return next_or_nothing(frame, pc)), pc3
        # Methods like f(x::Ref{<:Real}) that use gensymmed typevars will not have the *exact*
        # signature of the active method. So let's get the active signature.
        frame.pc = pc
        pc = define ? step_expr!(recurse, frame, stmt, true) : next_or_nothing!(frame)
        meth = whichtt(sigt)
        isa(meth, Method) && push!(signatures, meth.sig) # inner methods are not visible
        name == name3 && return pc, pc3     # if this was an inner method we should keep going
        stmt = pc_expr(frame, pc)  # there *should* be more statements in this frame
    end
end
methoddef!(@nospecialize(recurse), signatures, frame::Frame, pc::Int; define=true) =
    methoddef!(recurse, signatures, frame, pc_expr(frame, pc), pc; define=define)
function methoddef!(@nospecialize(recurse), signatures, frame::Frame; define=true)
    pc = frame.pc
    stmt = pc_expr(frame, pc)
    if !ismethod(stmt)
        pc = next_until!(ismethod, recurse, frame, true)
    end
    pc === nothing && error("pc at end of frame without finding a method")
    methoddef!(recurse, signatures, frame, pc; define=define)
end
methoddef!(signatures, frame::Frame; define=true) =
    methoddef!(finish_and_return!, signatures, frame; define=define)

function methoddefs!(@nospecialize(recurse), signatures, frame::Frame, pc; define=true)
    ret = methoddef!(recurse, signatures, frame, pc; define=define)
    pc = ret === nothing ? ret : ret[1]
    return _methoddefs!(recurse, signatures, frame, pc; define=define)
end
function methoddefs!(@nospecialize(recurse), signatures, frame::Frame; define=true)
    ret = methoddef!(recurse, signatures, frame; define=define)
    pc = ret === nothing ? ret : ret[1]
    return _methoddefs!(recurse, signatures, frame, pc; define=define)
end
methoddefs!(signatures, frame::Frame; define=true) =
    methoddefs!(finish_and_return!, signatures, frame; define=define)

function _methoddefs!(@nospecialize(recurse), signatures, frame::Frame, pc; define=define)
    while pc !== nothing
        stmt = pc_expr(frame, pc)
        if !ismethod(stmt)
            pc = next_until!(ismethod, recurse, frame, true)
        end
        pc === nothing && break
        ret = methoddef!(recurse, signatures, frame, pc; define=define)
        pc = ret === nothing ? ret : ret[1]
    end
    return pc
end

"""
    mbody = bodymethod(m::Method)

Return the "body method" for a method `m`. `mbody` contains the code of the function body
when `m` was defined.
"""
function bodymethod(mkw::Method)
    function is_self_call(stmt, slotnames, argno=1)
        if isa(stmt, Expr)
            if stmt.head == :call
                a = stmt.args[argno]
                if isa(a, SlotNumber)
                    if slotnames[a.id] == Symbol("#self#")
                        return true
                    end
                end
            end
        end
        return false
    end
    m = mkw
    local src
    while true
        framecode = JuliaInterpreter.get_framecode(m)
        fakeargs = Any[nothing for i = 1:length(framecode.scope.nargs)]
        frame = JuliaInterpreter.prepare_frame(framecode, fakeargs, isa(m.sig, UnionAll) ? sparam_ub(m) : Core.svec())
        src = framecode.src
        (length(src.code) > 1 && is_self_call(src.code[end-1], src.slotnames)) || break
        # Build the optional arg, so we can get its type
        pc = frame.pc
        while pc < length(src.code) - 1
            pc = step_expr!(frame)
        end
        val = pc > 1 ? frame.framedata.ssavalues[pc-1] : src.code[1].args[end]
        sig = Tuple{Base.unwrap_unionall(m.sig).parameters..., typeof(val)}
        m = whichtt(sig)
    end
    length(src.code) > 1 || return m
    stmt = src.code[end-1]
    if isexpr(stmt, :call) && (f = stmt.args[1]; isa(f, QuoteNode))
        # Check that it has a #self# call
        hasself = any(i->is_self_call(stmt, src.slotnames, i), 1:length(stmt.args))
        hasself || return m
        f = f.value
        mths = methods(f)
        if length(mths) == 1
            return first(mths)
        end
    end
    return m
end

function sparam_ub(meth::Method)
    typs = []
    sig = meth.sig
    while sig isa UnionAll
        push!(typs, Symbol(sig.var.ub))
        sig = sig.body
    end
    return Core.svec(typs...)
end

# precompilation

if ccall(:jl_generating_output, Cint, ()) == 1
    kwdefine = NamedTuple{(:define,),Tuple{Bool}}
    for ct in (Vector{Any}, Set{Any})
        f = methoddef!
        precompile(Tuple{typeof(f), ct, Any, Frame, Expr, Int})
        precompile(Tuple{Core.kwftype(typeof(f)), kwdefine, typeof(f), ct, Any, Frame, Expr, Int})
        f = methoddefs!
        precompile(Tuple{typeof(f), ct, Any, Frame})
        precompile(Tuple{Core.kwftype(typeof(f)), kwdefine, typeof(f), ct, Any, Frame})
    end
    precompile(Tuple{typeof(get_parentname), Symbol})
    precompile(Tuple{typeof(correct_name!), Any, Frame, Int, Symbol, Symbol})
end

end # module