text:
# abstract type Constraint end
"""
    Dirichlet(u, ∂Ω, f)
    Dirichlet(u, ∂Ω, f, component)

Create a Dirichlet boundary condition on `u` on the `∂Ω` part of
the boundary. `f` is a function that takes two arguments, `x` and `t`
where `x` is the spatial coordinate and `t` is the current time,
and returns the prescribed value. For example, here we create a
Dirichlet condition for the `:u` field, on the faceset called
`∂Ω` and the value given by the `sin` function:

```julia
dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t))
```

If `:u` is a vector field we can specify which component the condition
should be applied to by specifying `component`. `component` can be given
either as an integer, or as a vector, for example:

```julia
dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t), 1)      # applied to component 1
dbc = Dirichlet(:u, ∂Ω, (x, t) -> sin(t), [1, 3]) # applied to component 1 and 3
```

`Dirichlet` boundary conditions are added to a [`ConstraintHandler`](@ref)
which applies the condition via [`apply!`](@ref) and/or [`apply_zero!`](@ref).
"""
struct Dirichlet # <: Constraint
    f::Function # f(x,t) -> value
    faces::Union{Set{Int},Set{FaceIndex},Set{EdgeIndex},Set{VertexIndex}}
    field_name::Symbol
    components::Vector{Int} # components of the field
    local_face_dofs::Vector{Int}
    local_face_dofs_offset::Vector{Int}
end
function Dirichlet(field_name::Symbol, faces::Union{T}, f::Function, component::Int=1) where T
    Dirichlet(field_name, copy(faces), f, [component])
end
function Dirichlet(field_name::Symbol, faces::Set{T}, f::Function, components::AbstractVector{Int}) where T
    unique(components) == components || error("components not unique: $components")
    # issorted(components) || error("components not sorted: $components")
    return Dirichlet(f, copy(faces), field_name, Vector(components), Int[], Int[])
end

"""
    AffineConstraint(constrained_dofs::Int, master_dofs::Vector{Int}, coeffs::Vector{T}, b::T) where T

Define an affine/linear constraint to constrain dofs of the form `u_i = ∑(u[j] * a[j]) + b`.
"""
struct AffineConstraint{T}
    constrained_dof::Int
    entries::Vector{Pair{Int, T}} # masterdofs and factors
    b::T # inhomogeneity
end

"""
    ConstraintHandler

Collection of constraints.
"""
struct ConstraintHandler{DH<:AbstractDofHandler,T}
    dbcs::Vector{Dirichlet}
    acs::Vector{AffineConstraint}
    prescribed_dofs::Vector{Int}
    free_dofs::Vector{Int}
    inhomogeneities::Vector{T}
    dofmapping::Dict{Int,Int} # global dof -> index into dofs and inhomogeneities
    bcvalues::Vector{BCValues{T}}
    dh::DH
    closed::ScalarWrapper{Bool}
end

function ConstraintHandler(dh::AbstractDofHandler)
    @assert isclosed(dh)
    ConstraintHandler(Dirichlet[], AffineConstraint[], Int[], Int[], Float64[], Dict{Int,Int}(), BCValues{Float64}[], dh, ScalarWrapper(false))
end

"""
    RHSData

Stores the constrained columns and mean of the diagonal of stiffness matrix `A`.
"""
struct RHSData{T}
    m::T
    constrained_columns::SparseMatrixCSC{T, Int}
end

"""
    get_rhs_data(ch::ConstraintHandler, A::SparseMatrixCSC) -> RHSData

Returns the needed [`RHSData`](@ref) for [`apply_rhs!`](@ref).

This must be used when the same stiffness matrix is reused for multiple steps,
for example when timestepping, with different non-homogeneouos Dirichlet boundary
conditions.
"""
function get_rhs_data(ch::ConstraintHandler, A::SparseMatrixCSC)
    m = meandiag(A)
    constrained_columns = A[:, ch.prescribed_dofs]
    return RHSData(m, constrained_columns)
end

"""
    apply_rhs!(data::RHSData, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false)

Applies the boundary condition to the right-hand-side vector without modifying the stiffness matrix.

See also: [`get_rhs_data`](@ref).
"""
function apply_rhs!(data::RHSData, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false)
    K = data.constrained_columns
    @assert length(f) == size(K, 1)
    @boundscheck checkbounds(f, ch.prescribed_dofs)
    m = data.m

    # TODO: Can the loops be combined or does the order matter?
    @inbounds for i in 1:length(ch.inhomogeneities)
        v = ch.inhomogeneities[i]
        if !applyzero && v != 0
            for j in nzrange(K, i)
                f[K.rowval[j]] -= v * K.nzval[j]
            end
        end
    end
    @inbounds for ac in ch.acs
        global_dof = ac.constrained_dof
        for (d, v) in ac.entries
            f[d] += f[global_dof] * v
        end
        f[global_dof] = 0
    end
    @inbounds for i in 1:length(ch.inhomogeneities)
        d = ch.prescribed_dofs[i]
        v = ch.inhomogeneities[i]
        vz = applyzero ? zero(eltype(f)) : v
        f[d] = vz * m
    end
end

function Base.show(io::IO, ::MIME"text/plain", ch::ConstraintHandler)
    println(io, "ConstraintHandler:")
    if !isclosed(ch)
        print(io, "  Not closed!")
    else
        print(io, "  BCs:")
        for dbc in ch.dbcs
            print(io, "\n    ", "Field: ", dbc.field_name, ", ", "Components: ", dbc.components)
        end
    end
end

isclosed(ch::ConstraintHandler) = ch.closed[]
free_dofs(ch::ConstraintHandler) = ch.free_dofs
prescribed_dofs(ch::ConstraintHandler) = ch.prescribed_dofs

"""
    close!(ch::ConstraintHandler)

Close and finalize the `ConstraintHandler`.
"""
function close!(ch::ConstraintHandler)
    @assert(!isclosed(ch))

    # Make ch.prescribed_dofs unique and sorted, and do the same operations for ch.inhomogeneities
    # TODO: This is probably quite slow atm, and the unique!() and sort() functions can be combined?
    dofs_vals = unique(first, zip(ch.prescribed_dofs, ch.inhomogeneities))
    copy!(ch.prescribed_dofs, getindex.(dofs_vals, 1))
    copy!(ch.inhomogeneities, getindex.(dofs_vals, 2))

    I = sortperm(ch.prescribed_dofs)
    ch.prescribed_dofs .= ch.prescribed_dofs[I]
    ch.inhomogeneities .= ch.inhomogeneities[I]

    copy!(ch.free_dofs, setdiff(1:ndofs(ch.dh), ch.prescribed_dofs))

    for i in 1:length(ch.prescribed_dofs)
        ch.dofmapping[ch.prescribed_dofs[i]] = i
    end

    # TODO:
    # Do a bunch of checks to see if the affine constraints are linearly indepented etc.
    # If they are not, it is possible to automatically reformulate the constraints
    # such that they become independent. However, at this point, it is left to
    # the user to assure this.
    sort!(ch.acs, by = ac -> ac.constrained_dof)

    ch.closed[] = true
    return ch
end

function dbc_check(ch::ConstraintHandler, dbc::Dirichlet)
    # check input
    dbc.field_name in getfieldnames(ch.dh) || throw(ArgumentError("field $(dbc.field_name) does not exist in DofHandler, existing fields are $(getfieldnames(ch.dh))"))
    #TODO FIX!!
    #for component in dbc.components
    #    0 < component <= ndim(ch.dh, dbc.field_name) || error("component $component is not within the range of field $field which has $(ndim(ch.dh, field)) dimensions")
    #end
    if length(dbc.faces) == 0
        @warn("added Dirichlet Boundary Condition to set containing 0 entities")
    end
end

"""
    add!(ch::ConstraintHandler, dbc::Dirichlet)

Add a `Dirichlet` boundary condition to the `ConstraintHandler`.
"""
function add!(ch::ConstraintHandler, dbc::Dirichlet)
    dbc_check(ch, dbc)
    celltype = getcelltype(ch.dh.grid)
    @assert isconcretetype(celltype)

    field_idx = find_field(ch.dh, dbc.field_name)
    # Extract stuff for the field
    interpolation = getfieldinterpolation(ch.dh, field_idx)#ch.dh.field_interpolations[field_idx]
    field_dim = getfielddim(ch.dh, field_idx)#ch.dh.field_dims[field_idx] # TODO: I think we don't need to extract these here ...

    if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset
        bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument
    else
        bcvalue = BCValues(interpolation, default_interpolation(celltype), eltype(dbc.faces))
    end
    _add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(ch.dh, dbc.field_name), bcvalue)

    return ch
end

"""
    add!(ch::ConstraintHandler, ac::AffineConstraint)

Add the `AffineConstraint` to the `ConstraintHandler`.
"""
function add!(ch::ConstraintHandler, newac::AffineConstraint)
    # Basic error checking
    for ac in ch.acs
        (ac.constrained_dof == newac.constrained_dof) &&
            error("Constraint already exist for dof $(ac.constrained_dof)")
        any(x -> x.first == newac.constrained_dof, ac.entries) &&
            error("New constrained dof $(newac.constrained_dof) is already used as a master dof.")
    end
    push!(ch.acs, newac)
    push!(ch.prescribed_dofs, newac.constrained_dof)
    push!(ch.inhomogeneities, newac.b)
    return ch
end

function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcfaces::Set{Index}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid))) where {Index<:BoundaryIndex}
    local_face_dofs, local_face_dofs_offset =
        _local_face_dofs_for_bc(interpolation, field_dim, dbc.components, offset, boundaryfunction(eltype(bcfaces)))
    copy!(dbc.local_face_dofs, local_face_dofs)
    copy!(dbc.local_face_dofs_offset, local_face_dofs_offset)

    # loop over all the faces in the set and add the global dofs to `constrained_dofs`
    constrained_dofs = Int[]
    redundant_faces = Index[]
    for (cellidx, faceidx) in bcfaces
        if cellidx ∉ cellset
            push!(redundant_faces, Index(cellidx, faceidx)) # will be removed from dbc
            continue # skip faces that are not part of the cellset
        end
        _celldofs = fill(0, ndofs_per_cell(ch.dh, cellidx))
        celldofs!(_celldofs, ch.dh, cellidx) # extract the dofs for this cell
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        append!(constrained_dofs, _celldofs[local_face_dofs[r]]) # TODO: for-loop over r and simply push! to ch.prescribed_dofs
        @debug println("adding dofs $(_celldofs[local_face_dofs[r]]) to dbc")
    end

    setdiff!(dbc.faces, redundant_faces)
    # save it to the ConstraintHandler
    push!(ch.dbcs, dbc)
    push!(ch.bcvalues, bcvalue)
    append!(ch.prescribed_dofs, constrained_dofs)
    for _ in 1:length(constrained_dofs)
        push!(ch.inhomogeneities, NaN)
    end
    return ch
end

# Calculate which local dof index live on each face:
# face `i` have dofs `local_face_dofs[local_face_dofs_offset[i]:local_face_dofs_offset[i+1]-1]
function _local_face_dofs_for_bc(interpolation, field_dim, components, offset, boundaryfunc::F=faces) where F
    local_face_dofs = Int[]
    local_face_dofs_offset = Int[1]
    for (_, face) in enumerate(boundaryfunc(interpolation))
        for fdof in face, d in 1:field_dim
            if d in components
                push!(local_face_dofs, (fdof-1)*field_dim + d + offset)
            end
        end
        push!(local_face_dofs_offset, length(local_face_dofs) + 1)
    end
    return local_face_dofs, local_face_dofs_offset
end

function _add!(ch::ConstraintHandler, dbc::Dirichlet, bcnodes::Set{Int}, interpolation::Interpolation, field_dim::Int, offset::Int, bcvalue::BCValues, cellset::Set{Int}=Set{Int}(1:getncells(ch.dh.grid)))
    if interpolation !== default_interpolation(typeof(ch.dh.grid.cells[first(cellset)]))
        @warn("adding constraint to nodeset is not recommended for sub/super-parametric approximations.")
    end

    ncomps = length(dbc.components)
    nnodes = getnnodes(ch.dh.grid)
    interpol_points = getnbasefunctions(interpolation)
    _celldofs = fill(0, ndofs_per_cell(ch.dh, first(cellset)))
    node_dofs = zeros(Int, ncomps, nnodes)
    visited = falses(nnodes)
    for cell in CellIterator(ch.dh, collect(cellset)) # only go over cells that belong to current FieldHandler
        celldofs!(_celldofs, cell) # update the dofs for this cell
        for idx in 1:min(interpol_points, length(cell.nodes))
            node = cell.nodes[idx]
            if !visited[node]
                noderange = (offset + (idx-1)*field_dim + 1):(offset + idx*field_dim) # the dofs in this node
                for (i,c) in enumerate(dbc.components)
                    node_dofs[i,node] = _celldofs[noderange[c]]
                    @debug println("adding dof $(_celldofs[noderange[c]]) to node_dofs")
                end
                visited[node] = true
            end
        end
    end

    constrained_dofs = Int[]
    sizehint!(constrained_dofs, ncomps*length(bcnodes))
    sizehint!(dbc.local_face_dofs, length(bcnodes))
    for node in bcnodes
        if !visited[node]
            # either the node belongs to another field handler or it does not have dofs in the constrained field
            continue
        end
        for i in 1:ncomps
            push!(constrained_dofs, node_dofs[i,node])
        end
        push!(dbc.local_face_dofs, node) # use this field to store the node idx for each node
    end

    # save it to the ConstraintHandler
    copy!(dbc.local_face_dofs_offset, constrained_dofs) # use this field to store the global dofs
    push!(ch.dbcs, dbc)
    push!(ch.bcvalues, bcvalue)
    append!(ch.prescribed_dofs, constrained_dofs)
    for _ in 1:length(constrained_dofs)
        push!(ch.inhomogeneities, NaN)
    end
    return ch
end

# Updates the DBC's to the current time `time`
function update!(ch::ConstraintHandler, time::Real=0.0)
    @assert ch.closed[]
    for (i,dbc) in enumerate(ch.dbcs)
        # Function barrier
        _update!(ch.inhomogeneities, dbc.f, dbc.faces, dbc.field_name, dbc.local_face_dofs, dbc.local_face_dofs_offset,
                 dbc.components, ch.dh, ch.bcvalues[i], ch.dofmapping, convert(Float64, time))
    end
end

# for faces
function _update!(inhomogeneities::Vector{Float64}, f::Function, faces::Set{<:BoundaryIndex}, field::Symbol, local_face_dofs::Vector{Int}, local_face_dofs_offset::Vector{Int},
                  components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues,
                  dofmapping::Dict{Int,Int}, time::T) where {T}

    dim = getdim(dh.grid)
    _tmp_cellid = first(faces)[1]

    N = nnodes_per_cell(dh.grid, _tmp_cellid)
    xh = zeros(Vec{dim, T}, N) # pre-allocate
    _celldofs = fill(0, ndofs_per_cell(dh, _tmp_cellid))

    for (cellidx, faceidx) in faces
        cellcoords!(xh, dh, cellidx)
        celldofs!(_celldofs, dh, cellidx) # update global dofs for this cell

        # no need to reinit!, enough to update current_face since we only need geometric shape functions M
        facevalues.current_face[] = faceidx

        # local dof-range for this face
        r = local_face_dofs_offset[faceidx]:(local_face_dofs_offset[faceidx+1]-1)
        counter = 1

        for location in 1:getnquadpoints(facevalues)
            x = spatial_coordinate(facevalues, location, xh)
            bc_value = f(x, time)
            @assert length(bc_value) == length(components)

            for i in 1:length(components)
                # find the global dof
                globaldof = _celldofs[local_face_dofs[r[counter]]]
                counter += 1

                dbc_index = dofmapping[globaldof]
                inhomogeneities[dbc_index] = bc_value[i]
                @debug println("prescribing value $(bc_value[i]) on global dof $(globaldof)")
            end
        end
    end
end

# for nodes
function _update!(inhomogeneities::Vector{Float64}, f::Function, nodes::Set{Int}, field::Symbol, nodeidxs::Vector{Int}, globaldofs::Vector{Int},
                  components::Vector{Int}, dh::AbstractDofHandler, facevalues::BCValues,
                  dofmapping::Dict{Int,Int}, time::Float64)
    counter = 1
    for (idx, nodenumber) in enumerate(nodeidxs)
        x = dh.grid.nodes[nodenumber].x
        bc_value = f(x, time)
        @assert length(bc_value) == length(components)
        for v in bc_value
            globaldof = globaldofs[counter]
            counter += 1
            dbc_index = dofmapping[globaldof]
            inhomogeneities[dbc_index] = v
            @debug println("prescribing value $(v) on global dof $(globaldof)")
        end
    end
end

# Saves the dirichlet boundary conditions to a vtkfile.
# Values will have a 1 where bcs are active and 0 otherwise
function WriteVTK.vtk_point_data(vtkfile, ch::ConstraintHandler)
    unique_fields = []
    for dbc in ch.dbcs
        push!(unique_fields, dbc.field_name)
    end
    unique!(unique_fields)

    for field in unique_fields
        nd = ndim(ch.dh, field)
        data = zeros(Float64, nd, getnnodes(ch.dh.grid))
        for dbc in ch.dbcs
            dbc.field_name != field && continue
            if eltype(dbc.faces) <: BoundaryIndex
                functype = boundaryfunction(eltype(dbc.faces))
                for (cellidx, faceidx) in dbc.faces
                    for facenode in functype(ch.dh.grid.cells[cellidx])[faceidx]
                        for component in dbc.components
                            data[component, facenode] = 1
                        end
                    end
                end
            else
                for nodeidx in dbc.faces
                    for component in dbc.components
                        data[component, nodeidx] = 1
                    end
                end
            end
        end
        vtk_point_data(vtkfile, data, string(field, "_bc"))
    end
    return vtkfile
end

"""

    apply!(K::SparseMatrixCSC, rhs::AbstractVector, ch::ConstraintHandler)

Adjust the matrix `K` and right hand side `rhs` to account for the Dirichlet boundary
conditions specified in `ch` such that `K \\ rhs` gives the expected solution.

    apply!(v::AbstractVector, ch::ConstraintHandler)

Apply Dirichlet boundary conditions, specified in `ch`, to the solution vector `v`.

# Examples
```julia
K, f = assemble_system(...) # Assemble system
apply!(K, f, ch)            # Adjust K and f to account for boundary conditions
u = K \\ f                   # Solve the system, u should be "approximately correct"
apply!(u, ch)               # Explicitly make sure bcs are correct
```

!!! note
    The last operation is not strictly necessary since the boundary conditions should
    already be fulfilled after `apply!(K, f, ch)`. However, solvers of linear systems are
    not exact, and thus `apply!(u, ch)` can be used to make sure the boundary conditions
    are fulfilled exactly.
"""
apply!

"""
    apply_zero!(K::SparseMatrixCSC, rhs::AbstractVector, ch::ConstraintHandler)

Adjust the matrix `K` and the right hand side `rhs` to account for prescribed Dirichlet
boundary conditions such that `du = K \\ rhs` give the expected result (e.g. with `du` zero
for all prescribed degrees of freedom).

    apply_zero!(v::AbstractVector, ch::ConstraintHandler)

Zero-out values in `v` corresponding to prescribed degrees of freedom.

These methods are typically used in e.g. a Newton solver where the increment, `du`, should
be prescribed to zero even for non-homogeneouos boundary conditions.

See also: [`apply!`](@ref).

# Examples
```julia
u = un + Δu                 # Current guess
K, g = assemble_system(...) # Assemble residual and tangent for current guess
apply_zero!(K, g, ch)       # Adjust tangent and residual to take prescribed values into account
ΔΔu = - K \\ g               # Compute the increment, prescribed values are "approximately" zero
apply_zero!(ΔΔu, ch)        # Make sure values are exactly zero
Δu .+= ΔΔu                  # Update current guess
```

!!! note
    The last call to `apply_zero!` is not strictly necessary since the boundary conditions
    should already be fulfilled after `apply!(K, g, ch)`. However, solvers of linear
    systems are not exact, and thus `apply!(ΔΔu, ch)` can be used to make sure the values
    for the prescribed degrees of freedom are fulfilled exactly.
"""
apply_zero!

apply_zero!(v::AbstractVector, ch::ConstraintHandler) = _apply_v(v, ch, true)
apply!(     v::AbstractVector, ch::ConstraintHandler) = _apply_v(v, ch, false)

function _apply_v(v::AbstractVector, ch::ConstraintHandler, apply_zero::Bool)
    @assert length(v) >= ndofs(ch.dh)
    v[ch.prescribed_dofs] .= apply_zero ? 0.0 : ch.inhomogeneities
    # Apply affine constraints, e.g u2 = u6 + b
    for ac in ch.acs
        for (d, s) in ac.entries
            v[ac.constrained_dof] += s * v[d]
        end
    end
    return v
end

function apply!(K::Union{SparseMatrixCSC,Symmetric}, ch::ConstraintHandler)
    apply!(K, eltype(K)[], ch, true)
end

function apply_zero!(K::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler)
    apply!(K, f, ch, true)
end

@enumx ApplyStrategy Transpose Inplace
# For backwards compatibility
const APPLY_TRANSPOSE = ApplyStrategy.Transpose
const APPLY_INPLACE = ApplyStrategy.Inplace

function apply!(KK::Union{SparseMatrixCSC,Symmetric}, f::AbstractVector, ch::ConstraintHandler, applyzero::Bool=false;
                strategy::ApplyStrategy.T=ApplyStrategy.Transpose)
    K = isa(KK, Symmetric) ? KK.data : KK
    @assert length(f) == 0 || length(f) == size(K, 1)
    @boundscheck checkbounds(K, ch.prescribed_dofs, ch.prescribed_dofs)
    @boundscheck length(f) == 0 || checkbounds(f, ch.prescribed_dofs)

    m = meandiag(K) # Use the mean of the diagonal here to not ruin things for iterative solver

    # Add inhomogeneities to f: (f - K * ch.inhomogeneities)
    if !applyzero
        @inbounds for i in 1:length(ch.inhomogeneities)
            d = ch.prescribed_dofs[i]
            v = ch.inhomogeneities[i]
            if v != 0
                for j in nzrange(K, d)
                    f[K.rowval[j]] -= v * K.nzval[j]
                end
            end
        end
    end

    # Condense K (C' * K * C) and f (C' * f)
    _condense!(K, f, ch.acs)

    # Remove constrained dofs from the matrix
    zero_out_columns!(K, ch.prescribed_dofs)
    if strategy === ApplyStrategy.Transpose
        K′ = copy(K)
        transpose!(K′, K)
        zero_out_columns!(K′, ch.prescribed_dofs)
        transpose!(K, K′)
    elseif strategy === ApplyStrategy.Inplace
        K[ch.prescribed_dofs, :] .= 0
    else
        error("Unknown apply strategy")
    end

    # Add meandiag to constraint dofs
    @inbounds for i in 1:length(ch.inhomogeneities)
        d = ch.prescribed_dofs[i]
        K[d, d] = m
        if length(f) != 0
            vz = applyzero ? zero(eltype(f)) : ch.inhomogeneities[i]
            f[d] = vz * m
        end
    end
end

# Similar to Ferrite._condense!(K, ch), but only add the non-zero entries to K (that arises from the condensation process)
function _condense_sparsity_pattern!(K::SparseMatrixCSC, acs::Vector{AffineConstraint})
    ndofs = size(K, 1)

    # Store linear constraint index for each constrained dof
    distribute = Dict{Int,Int}(acs[c].constrained_dof => c for c in 1:length(acs))

    for col in 1:ndofs
        # Since we will possibly be pushing new entries to K, the field K.rowval will grow.
        # Therefor we must extract this before iterating over K
        range = nzrange(K, col)
        _rows = K.rowval[range]
        dcol = get(distribute, col, 0)
        if dcol == 0
            for row in _rows
                drow = get(distribute, row, 0)
                if drow != 0
                    ac = acs[drow]
                    for (d, _) in ac.entries
                        add_entry!(K, d, col)
                    end
                end
            end
        else
            for row in _rows
                drow = get(distribute, row, 0)
                if drow == 0
                    ac = acs[dcol]
                    for (d, _) in ac.entries
                        add_entry!(K, row, d)
                    end
                else
                    ac1 = acs[dcol]
                    for (d1, _) in ac1.entries
                        ac2 = acs[distribute[row]]
                        for (d2, _) in ac2.entries
                            add_entry!(K, d1, d2)
                        end
                    end
                end
            end
        end
    end
end

# Condenses K and f: C'*K*C, C'*f, in-place assuming the sparsity pattern is correct
function _condense!(K::SparseMatrixCSC, f::AbstractVector, acs::Vector{AffineConstraint})

    ndofs = size(K, 1)
    condense_f = !(length(f) == 0)
    condense_f && @assert( length(f) == ndofs )

    # Store linear constraint index for each constrained dof
    distribute = Dict{Int,Int}(acs[c].constrained_dof => c for c in 1:length(acs))

    for col in 1:ndofs
        dcol = get(distribute, col, 0)
        if dcol == 0
            for a in nzrange(K, col)
                row = K.rowval[a]
                drow = get(distribute, row, 0)
                if drow != 0
                    ac = acs[drow]
                    for (d, v) in ac.entries
                        Kval = K.nzval[a]
                        _addindex_sparsematrix!(K, v * Kval, d, col)
                    end

                    # Perform f - K*g. However, this has already been done in outside this functions so we skip this.
                    #if condense_f
                    #    f[col] -= K.nzval[a] * ac.b;
                    #end
                end
            end
        else
            for a in nzrange(K, col)
                row = K.rowval[a]
                drow = get(distribute, row, 0)
                if drow == 0
                    ac = acs[dcol]
                    for (d,v) in ac.entries
                        Kval = K.nzval[a]
                        _addindex_sparsematrix!(K, v * Kval, row, d)
                    end
                else
                    ac1 = acs[dcol]
                    for (d1,v1) in ac1.entries
                        ac2 = acs[drow]
                        for (d2,v2) in ac2.entries
                            Kval = K.nzval[a]
                            _addindex_sparsematrix!(K, v1 * v2 * Kval, d1, d2)
                        end
                    end
                end
            end

            if condense_f
                ac = acs[dcol]
                for (d,v) in ac.entries
                    f[d] += f[col] * v
                end
                @assert ac.constrained_dof == col
                f[ac.constrained_dof] = 0.0
            end
        end
    end
end

# Copied from SparseArrays._setindex_scalar!(...)
# Custom SparseArrays._setindex_scalar!() that throws error if entry K(_i,_j) does not exist
function _addindex_sparsematrix!(A::SparseMatrixCSC{Tv,Ti}, v::Tv, i::Ti, j::Ti) where {Tv, Ti}
    if !((1 <= i <= size(A, 1)) & (1 <= j <= size(A, 2)))
        throw(BoundsError(A, (i,j)))
    end
    coljfirstk = Int(SparseArrays.getcolptr(A)[j])
    coljlastk = Int(SparseArrays.getcolptr(A)[j+1] - 1)
    searchk = searchsortedfirst(rowvals(A), i, coljfirstk, coljlastk, Base.Order.Forward)
    if searchk <= coljlastk && rowvals(A)[searchk] == i
        # Column j contains entry A[i,j]. Update and return
        nonzeros(A)[searchk] += v
        return A
    end
    error("Sparsity pattern missing entries for the condensation pattern. Make sure to call `create_sparsity_pattern(dh::DofHandler, ch::ConstraintHandler) when using linear constraints.`")
end

# A[i,j] += 0.0 does not add entries to sparse matrices, so we need to first add 1.0, and then remove it
# TODO: Maybe this can be done for vectors i and j instead of doing it individually?
function add_entry!(A::SparseMatrixCSC, i::Int, j::Int)
    if iszero(A[i,j]) # Check first if zero to not remove already non-zero entries
        A[i,j] = 1
        A[i,j] = 0
    end
end

"""
    create_constraint_matrix(ch::ConstraintHandler)

Create and return the constraint matrix, `C`, and the inhomogeneities, `g`, from the affine
(linear) and Dirichlet constraints in `ch`.

The constraint matrix relates constrained, `a_c`, and free, `a_f`, degrees of freedom via
`a_c = C * a_f + g`. The condensed system of linear equations is obtained as
`C' * K * C = C' *  (f - K * g)`.
"""
function create_constraint_matrix(ch::ConstraintHandler{dh,T}) where {dh,T}
    @assert(isclosed(ch))

    I = Int[]; J = Int[]; V = T[];
    g = zeros(T, ndofs(ch.dh)) # inhomogeneities

    for (j, d) in enumerate(ch.free_dofs)
       push!(I, d)
       push!(J, j)
       push!(V, 1.0)
    end

    for ac in ch.acs
        for (d, v) in ac.entries
            push!(I, ac.constrained_dof)
            j = searchsortedfirst(ch.free_dofs, d)
            push!(J, j)
            push!(V, v)
        end
        g[ac.constrained_dof] = ac.b
    end
    g[ch.prescribed_dofs] .= ch.inhomogeneities

    C = sparse(I, J,  V, ndofs(ch.dh), length(ch.free_dofs))

    return C, g
end

# columns need to be stored entries, this is not checked
function zero_out_columns!(K, dofs::Vector{Int}) # can be removed in 0.7 with #24711 merged
    @debug @assert issorted(dofs)
    for col in dofs
        r = nzrange(K, col)
        K.nzval[r] .= 0.0
    end
end

function meandiag(K::AbstractMatrix)
    z = zero(eltype(K))
    for i in 1:size(K, 1)
        z += abs(K[i, i])
    end
    return z / size(K, 1)
end

#Function for adding constraint when using multiple celltypes
function add!(ch::ConstraintHandler, fh::FieldHandler, dbc::Dirichlet)
    _check_cellset_dirichlet(ch.dh.grid, fh.cellset, dbc.faces)

    celltype = getcelltype(ch.dh.grid, first(fh.cellset)) #Assume same celltype of all cells in fh.cellset

    field_idx = find_field(fh, dbc.field_name)
    # Extract stuff for the field
    interpolation = getfieldinterpolations(fh)[field_idx]
    field_dim = getfielddims(fh)[field_idx]

    if eltype(dbc.faces)==Int #Special case when dbc.faces is a nodeset
        bcvalue = BCValues(interpolation, default_interpolation(celltype), FaceIndex) #Not used by node bcs, but still have to pass it as an argument
    else
        bcvalue = BCValues(interpolation, default_interpolation(celltype), eltype(dbc.faces))
    end

    Ferrite._add!(ch, dbc, dbc.faces, interpolation, field_dim, field_offset(fh, dbc.field_name), bcvalue, fh.cellset)
    return ch
end

function _check_cellset_dirichlet(::AbstractGrid, cellset::Set{Int}, faceset::Set{<:BoundaryIndex})
    for (cellid,faceid) in faceset
        if !(cellid in cellset)
            @warn("You are trying to add a constraint to a face that is not in the cellset of the fieldhandler. The face will be skipped.")
        end
    end
end

function _check_cellset_dirichlet(grid::AbstractGrid, cellset::Set{Int}, nodeset::Set{Int})
    nodes = Set{Int}()
    for cellid in cellset
        for nodeid in grid.cells[cellid].nodes
            nodeid ∈ nodes || push!(nodes, nodeid)
        end
    end

    for nodeid in nodeset
        if !(nodeid ∈ nodes)
            @warn("You are trying to add a constraint to a node that is not in the cellset of the fieldhandler. The node will be skipped.")
        end
    end
end

create_symmetric_sparsity_pattern(dh::MixedDofHandler, ch::ConstraintHandler) = Symmetric(_create_sparsity_pattern(dh, ch, true), :U)
create_symmetric_sparsity_pattern(dh::DofHandler,      ch::ConstraintHandler) = Symmetric(_create_sparsity_pattern(dh, ch, true), :U)

create_sparsity_pattern(dh::MixedDofHandler, ch::ConstraintHandler) = _create_sparsity_pattern(dh, ch, false)
create_sparsity_pattern(dh::DofHandler,      ch::ConstraintHandler) = _create_sparsity_pattern(dh, ch, false)

"""
    PeriodicDirichlet(u, Γ⁻ => Γ⁺, component=1)
    PeriodicDirichlet(u, Γ⁻ => Γ⁺, f, component=1)

Create a periodic Dirichlet boundary condition for the field `u`, with a mirror boundary,
`Γ⁻` and an image boundary, `Γ⁺`. The condition is imposed in a strong sense, and requires
(i) a periodic domain (usually a cube) and (ii) a periodic mesh.

See also manual section on [Periodic boundary conditions](@ref).
"""
struct PeriodicDirichlet
    field_name::Symbol
    components::Vector{Int} # components of the field
    face_pairs::Vector{Pair{String,String}}
    func::Union{Function,Nothing}
end

PeriodicDirichlet(fn::Symbol, fp::Vector{<:Pair}, c::Union{Int,Vector{Int}}=1) =
    PeriodicDirichlet(fn, fp, nothing, c)

PeriodicDirichlet(fn::Symbol, fp::Vector{<:Pair}, f::Union{Function,Nothing}, c::Union{Int,Vector{Int}}=1) =
    PeriodicDirichlet(fn, sort!(vec(collect(c))), fp, f)

function add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet)
    field_idx = find_field(ch.dh, pdbc.field_name)
    interpolation = getfieldinterpolation(ch.dh, field_idx)
    field_dim = getfielddim(ch.dh, field_idx)
    _add!(ch, pdbc, interpolation, field_dim, field_offset(ch.dh, pdbc.field_name))
    return ch
end

function _add!(ch::ConstraintHandler, pdbc::PeriodicDirichlet, interpolation::Interpolation,
               field_dim::Int, offset::Int)
    grid = ch.dh.grid
    Tx = typeof(first(ch.dh.grid.nodes).x) # Vec{D,T}
    # TODO: This requires full periodicity for now. (Why? I don't remember,
    #       but something in the code below assumes this... I think...)
    @assert length(pdbc.face_pairs) == length(Tx)

    # Indices of the local dofs for the faces
    local_face_dofs, local_face_dofs_offset =
        _local_face_dofs_for_bc(interpolation, field_dim, pdbc.components, offset)
    mirrored_local_face_dofs =
        mirror_local_dofs(local_face_dofs, local_face_dofs_offset, interpolation, length(pdbc.components))

    # First, create the facemap which maps mirror faces to image faces
    face_map = Dict{FaceIndex,FaceIndex}() # mirror face => image face
    dof_map = Dict{Int,Int}() # mirror dof => image => dof

    # Some stuff for finding the corner points
    all_node_idxs = Set{Int}()
    max_x = Tx(i -> typemin(eltype(Tx)))
    min_x = Tx(i -> typemax(eltype(Tx)))

    for (xyz, (mirror, image)) in enumerate(pdbc.face_pairs)
        # Need the order here, so collect
        mirror_faceset = collect(getfaceset(grid, mirror))
        mirror_mean_x = Tx[] # mean face coordinates
        for (c, f) in mirror_faceset
            fn = faces(grid.cells[c])[f]
            push!(mirror_mean_x, sum(grid.nodes[i].x for i in fn) / length(fn))
            # Also keep track of all nodes for corner finding
            union!(all_node_idxs, grid.cells[c].nodes)
        end
        # Same dance for the image
        image_faceset = collect(getfaceset(grid,  image))
        image_mean_x = Tx[]
        for (c, f) in image_faceset
            fn = faces(grid.cells[c])[f]
            push!(image_mean_x, sum(grid.nodes[i].x for i in fn) / length(fn))
            # Also keep track of all nodes for corner finding
            union!(all_node_idxs, grid.cells[c].nodes)
        end
        function extract_single_x_on_face(fs) # :)
            c, f = first(fs)
            fn = faces(grid.cells[c])[f]
            return grid.nodes[fn[1]].x
        end
        # Add a shift to help the tree search
        xmirror = extract_single_x_on_face(mirror_faceset)
        ximage  = extract_single_x_on_face(image_faceset)
        xdist = ximage[xyz] - xmirror[xyz]
        @assert xdist > 0
        shift = Tx(i -> i == xyz ? xdist : 0)
        for i in eachindex(mirror_mean_x)
            mirror_mean_x[i] += shift
        end
        # Use KDTree to find closest face
        tree = KDTree(image_mean_x)
        idxs, _ = NearestNeighbors.nn(tree, mirror_mean_x)
        for (m, idx) in zip(mirror_faceset, idxs)
            face_map[m] = image_faceset[idx]
        end
        # Update maximum/minimum coords from the bounding box of the tree
        max_x = Tx(i -> i == xyz ? max(max_x[i], tree.hyper_rec.maxes[i]) : max_x[i])
        min_x = Tx(i -> i == xyz ? min(min_x[i], tree.hyper_rec.maxes[i] - xdist) : min_x[i])

        # Create the constraints
        mirror_dofs = zeros(Int, ndofs_per_cell(ch.dh))
         image_dofs = zeros(Int, ndofs_per_cell(ch.dh))
        for (m, i) in face_map
            celldofs!(mirror_dofs, ch.dh, m[1])
            celldofs!( image_dofs, ch.dh, i[1])

            mdof_range = local_face_dofs_offset[m[2]] : (local_face_dofs_offset[m[2] + 1] - 1)
            idof_range = local_face_dofs_offset[i[2]] : (local_face_dofs_offset[i[2] + 1] - 1)

            for (md, id) in zip(mdof_range, idof_range)
                cdof = mirror_dofs[mirrored_local_face_dofs[md]]
                mdof = image_dofs[local_face_dofs[id]]

                if haskey(dof_map, mdof)
                    mdof′ = dof_map[mdof]
                    # @info "$cdof => $mdof, but $mdof => $mdof′, remapping $cdof => $mdof′."
                    push!(dof_map, cdof => mdof′)
                # elseif haskey(dof_map, cdof) && dof_map[cdof] == mdof
                    # @info "$cdof => $mdof already in the set, skipping."
                elseif haskey(dof_map, cdof)
                    # @info "$cdof => $mdof, but $cdof => $(dof_map[cdof]) already, skipping."
                else
                    # @info "$cdof => $mdof."
                    push!(dof_map, cdof => mdof)
                end
            end
        end
    end

    # Find the corners and put homogeneous Dirichlet there instead
    all_node_idxs_v = collect(all_node_idxs)
    points = construct_cornerish(min_x, max_x)
    tree = KDTree(Tx[grid.nodes[i].x for i in all_node_idxs_v])
    idxs, _ = NearestNeighbors.nn(tree, points)
    corner_set = Set{Int}(all_node_idxs_v[i] for i in idxs)

    dbc = Dirichlet(pdbc.field_name, corner_set,
        pdbc.func === nothing ? (x, _) -> pdbc.components * eltype(x)(0) : pdbc.func,
        pdbc.components
    )

    # Create a temp constraint handler just to find the dofs in the nodes...
    chtmp = ConstraintHandler(ch.dh)
    add!(chtmp, dbc)
    close!(chtmp)
    # No need to update!(chtmp, t) here since we only care about the dofs
    # TODO: Right? maybe if the user passed f we need to...
    foreach(x -> delete!(dof_map, x), chtmp.prescribed_dofs)

    # Need to reset the internal of this DBC in order to add! it again...
    resize!(dbc.local_face_dofs, 0)
    resize!(dbc.local_face_dofs_offset, 0)

    # Create another temp constraint handler if we need to compute inhomogeneities
    inhomogeneity_map = nothing
    if pdbc.func !== nothing
        chtmp2 = ConstraintHandler(ch.dh)
        all_faces = union!(Set{FaceIndex}(),
            (getfaceset(grid, x) for x in Iterators.flatten(pdbc.face_pairs))...
        )
        dbc_all = Dirichlet(pdbc.field_name, all_faces, pdbc.func, pdbc.components)
        add!(chtmp2, dbc_all); close!(chtmp2)
        # Call update! here since we need it to construct the affine constraints...
        # TODO: This doesn't allow for time dependent constraints...
        update!(chtmp2, 0.0)
        inhomogeneity_map = Dict{Int, Float64}()
        for (k, v) in dof_map
            g = chtmp2.inhomogeneities
            push!(inhomogeneity_map,
                  v => - g[chtmp2.dofmapping[v]] + g[chtmp2.dofmapping[k]]
            )
        end
    end

    # Any remaining mappings are added as homogeneous AffineConstraints
    for (k, v) in dof_map
        ac = AffineConstraint(k, [v => 1.0], inhomogeneity_map === nothing ? 0.0 : inhomogeneity_map[v])
        add!(ch, ac)
    end
    # Add the Dirichlet for the corners
    add!(ch, dbc)

    return ch
end
function construct_cornerish(min_x::V, max_x::V) where {T, V <: Vec{1,T}}
    lx = max_x - min_x
    max_x += lx
    min_x -= lx
    return V[min_x, max_x]
end
function construct_cornerish(min_x::V, max_x::V) where {T, V <: Vec{2,T}}
    lx = max_x - min_x
    max_x += lx
    min_x -= lx
    return V[
       max_x,
       min_x,
       Vec{2,T}((max_x[1], min_x[2])),
       Vec{2,T}((min_x[1], max_x[2])),
    ]
end
function construct_cornerish(min_x::V, max_x::V) where {T, V <: Vec{3,T}}
    error("not implemented yet, please contribute :)")
end

function mirror_local_dofs(_, _, ::Lagrange{1}, ::Int)
    # For 1D there is nothing to do
end
function mirror_local_dofs(local_face_dofs, local_face_dofs_offset, ip::Lagrange{2,<:Union{RefCube,RefTetrahedron}}, n::Int)
    # For 2D we always permute since Ferrite defines dofs counter-clockwise
    ret = copy(local_face_dofs)
    for (i, f) in enumerate(faces(ip))
        this_offset = local_face_dofs_offset[i]
        other_offset = this_offset + n
        for d in 1:n
            idx1 = this_offset + (d - 1)
            idx2 = other_offset + (d - 1)
            tmp = ret[idx1]
            ret[idx1] = ret[idx2]
            ret[idx2] = tmp
        end
    end
    return ret
end
function mirror_local_dofs(_, _, ::Lagrange{3})
    error("not implemented yet, please contribute :)")
end


text:
module Linearisers

using ForwardDiff

export get_jacobian

"""
    get_jacobian(func, args...)

Compute the jacobian of given function `func`.
`args...` should be arguments of `func`.

# Equations
jacobian = [[∂f1/∂x1 ... ∂f1/∂xn];
                ⋮           ⋮
            [∂fm/∂x1 ... ∂fm/∂xn]]

# Examples
```julia
using SuccessiveConvexProgrammings
using LinearAlgebra

function my_func_array(x, u, t, k)
    A = [1 2;
         3 4;
         5 6]
    B = Matrix(I, 3, 3)
    return A*x + B*u
end

x = [1, 2]
u = [3, 4, 5]
t = [6, 7]
k = [8, 9]
jacob = get_jacobian(my_func_array, x, u, t, k)
# result
3×9 Array{Float64,2}:
 1.0  2.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
 3.0  4.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
 5.0  6.0  0.0  0.0  1.0  0.0  0.0  0.0  0.0
```
"""
function get_jacobian(func, args...)
    jacob = zeros(length(func(args...)), 0)  # dummy calculation
    args_dim = 0
    for (i, arg) in enumerate(args)
        args_tmp = [args...]
        args_tmp[1], args_tmp[i] = args_tmp[i], args_tmp[1]
        altered_func = function(args2...)
            args_tmp2 = [args2...]
            args_tmp2[1], args_tmp2[i] = args_tmp2[i], args_tmp2[1]
            return func(args_tmp2...)
        end
        new_func = function(arg)
            return altered_func(arg, args_tmp[2:end]...)
        end
        jacob = cat(jacob, ForwardDiff.jacobian(new_func, arg), dims=2)
    end
    return jacob
end

end

ext: jl


text:
## out-of-place operator backends
_GetDeriv(ADmode::Val{:ReverseDiff}; kwargs...) = throw("GetDeriv() not available for ReverseDiff.jl")
_GetGrad(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.gradient
_GetJac(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.jacobian
_GetHess(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.hessian

_GetDeriv(ADmode::Val{:Zygote}; kwargs...) = throw("GetDeriv() not available for Zygote.jl")
_GetGrad(ADmode::Val{:Zygote}; order::Int=-1, kwargs...) = (Func::Function,p;Kwargs...) -> Zygote.gradient(Func, p; kwargs...)[1]
_GetJac(ADmode::Val{:Zygote}; order::Int=-1, kwargs...) = (Func::Function,p;Kwargs...) -> Zygote.jacobian(Func, p; kwargs...)[1]
_GetHess(ADmode::Val{:Zygote}; order::Int=-1, kwargs...) = (Func::Function,p;Kwargs...) -> Zygote.hessian(Func, p; kwargs...)
_GetDoubleJac(ADmode::Val{:Zygote}; kwargs...) = throw("GetDoubleJac() not available for Zygote.jl") # Zygote does not support mutating arrays

_GetDeriv(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_derivative
_GetGrad(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_gradient
_GetJac(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_jacobian
_GetHess(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_hessian
_GetDoubleJac(ADmode::Val{:FiniteDiff}; kwargs...) = throw("GetDoubleJac() not available for FiniteDiff.jl")

_GetDeriv(ADmode::Val{:FiniteDifferences}; kwargs...) = throw("GetDeriv() not available for FiniteDifferences.jl")
_GetGrad(ADmode::Val{:FiniteDifferences}; order::Int=3, kwargs...) = (Func::Function,p;Kwargs...) -> FiniteDifferences.grad(central_fdm(order,1), Func, p; kwargs...)[1]
_GetJac(ADmode::Val{:FiniteDifferences}; order::Int=3, kwargs...) = (Func::Function,p;Kwargs...) -> FiniteDifferences.jacobian(central_fdm(order,1), Func, p; kwargs...)[1]
_GetHess(ADmode::Val{:FiniteDifferences}; order::Int=5, kwargs...) = (Func::Function,p;Kwargs...) -> FiniteDifferences.jacobian(central_fdm(order,1), z->FiniteDifferences.grad(central_fdm(order,1), Func, z)[1], p)[1]

## in-place operator backends
_GetGrad!(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.gradient!
_GetJac!(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.jacobian!
_GetHess!(ADmode::Val{:ReverseDiff}; kwargs...) = ReverseDiff.hessian!
_GetMatrixJac!(ADmode::Val{:ReverseDiff}; kwargs...) = _GetJac!(ADmode; kwargs...) # DELIBERATE!!!! _GetJac!() recognizes output format from given Array

#_GetDeriv!(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_derivative!
_GetGrad!(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_gradient!
function _GetJac!(ADmode::Val{:FiniteDiff}; kwargs...)
    function FiniteDiff__finite_difference_jacobian!(Y::AbstractArray{<:Number}, F::Function, X, args...; kwargs...)
        # in-place FiniteDiff operators assume that function itself is also in-place
        if MaximalNumberOfArguments(F) > 1
            FiniteDiff.finite_difference_jacobian!(Y, F, X, args...; kwargs...)
        else
            # Use fake method
            (Y[:] .= vec(_GetJac(ADmode; kwargs...)(F, X, args...)))
            # FiniteDiff.finite_difference_jacobian!(Y, (Res,x)->copyto!(Res,F(x)), args...; kwargs...)
        end
    end
end
_GetHess!(ADmode::Val{:FiniteDiff}; kwargs...) = FiniteDiff.finite_difference_hessian!
_GetMatrixJac!(ADmode::Val{:FiniteDiff}; kwargs...) = _GetJac!(ADmode; kwargs...)

# Fake in-place methods
function _GetGrad!(ADmode::Union{<:Val{:Zygote},<:Val{:FiniteDifferences}}; verbose::Bool=false, kwargs...)
    verbose && (@warn "Using fake in-place differentiation operator GetGrad!() for ADmode=$ADmode because backend does not supply appropriate method.")
    FakeInPlaceGrad!(Y::AbstractVector,F::Function,X::AbstractVector) = copyto!(Y, _GetGrad(ADmode; kwargs...)(F, X))
end
function _GetJac!(ADmode::Union{Val{:Zygote},<:Val{:FiniteDifferences}}; verbose::Bool=false, kwargs...)
    verbose && (@warn "Using fake in-place differentiation operator GetJac!() for ADmode=$ADmode because backend does not supply appropriate method.")
    FakeInPlaceJac!(Y::AbstractMatrix,F::Function,X::AbstractVector) = copyto!(Y, _GetJac(ADmode; kwargs...)(F, X))
end
function _GetHess!(ADmode::Union{Val{:Zygote},<:Val{:FiniteDifferences}}; verbose::Bool=false, kwargs...)
    verbose && (@warn "Using fake in-place differentiation operator GetHess!() for ADmode=$ADmode because backend does not supply appropriate method.")
    FakeInPlaceHess!(Y::AbstractMatrix,F::Function,X::AbstractVector) = copyto!(Y, _GetHess(ADmode; kwargs...)(F, X))
end
function _GetMatrixJac!(ADmode::Union{Val{:Zygote},<:Val{:FiniteDifferences}}; verbose::Bool=false, kwargs...)
    verbose && (@warn "Using fake in-place differentiation operator GetMatrixJac!() for ADmode=$ADmode because backend does not supply appropriate method.")
    FakeInPlaceMatrixJac!(Y::AbstractArray,F::Function,X::AbstractVector) = (Y[:] .= vec(_GetJac(ADmode; kwargs...)(F, X)))
end


text:
### Kernel Density Estimation
export
    kde,
    KDEMethod,
    FixedBandwidth,
    SilvermansRule,
    bandwidth
include("kde_methods.jl")

export
    gmm_matrix,
    gmm_matrix!
include("gmm_matrix.jl")

### Information Measure calculation
export
    information,
    InformationMeasure,
    NoInformation
export
    ParticleInformationMeasure,
    DiscreteEntropy,
    discrete_entropy
export
    KDEInformationMeasure,
    KDEEntropy,
    GMMEntropy

abstract type InformationMeasure end
"""
    Measures how much information is contained in the given belief.
    Measuring is done with the given InformationMeasure.
"""
information(m::InformationMeasure, args...; kwargs...)::Float64 = (error("information not implemented for $m"); return 0.0)

struct NoInformation <: InformationMeasure end
information(::NoInformation; kwargs...)::Float64 = 0.0

abstract type ParticleInformationMeasure <: InformationMeasure end

struct DiscreteEntropy <: ParticleInformationMeasure end
information(::DiscreteEntropy; bp::AbstractParticleBelief, kwargs...)::Float64 = -discrete_entropy(bp)
# Works only if enough particles in every discrete state, by summing up the weights
# for continuous distributions where all particles are in different states this yields -log(1/n_particles(b))
discrete_entropy(b::AbstractParticleBelief)::Float64 = -sum(p * log(p) for p in values(ParticleFilters.probdict(b)) if p > 0.0)

abstract type KDEInformationMeasure <: InformationMeasure end

struct GMMEntropy <: KDEInformationMeasure
    kde_method::KDEMethod
end
function information(m::GMMEntropy; bp::AbstractParticleBelief, kwargs...)::Float64
    ws = weight_sum(bp)
    dot(weights(bp)/ws, log.(gmm_matrix(bp, bandwidth(m.kde_method, bp)) * weights(bp)/ws))
end


text:
using BinaryBuilder

name = "SCS"
version = v"2.1.2"

# Collection of sources required to build SCSBuilder
sources = [
    GitSource("https://github.com/cvxgrp/scs.git", "4ed6c2abf28399c01a0417ff3456b2639560afa6")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/scs*
flags="DLONG=1 USE_OPENMP=0"
blasldflags="-L${prefix}/lib"
# see https://github.com/JuliaPackaging/Yggdrasil/blob/0bc1abd56fa176e3d2cc2e48e7bf85a26c948c40/OpenBLAS/build_tarballs.jl#L23
if [[ ${nbits} == 64 ]] && [[ ${target} != aarch64* ]]; then
    flags="${flags} BLAS64=1 BLASSUFFIX=_64_"
    blasldflags+=" -lopenblas64_"
else
    blasldflags+=" -lopenblas"
fi

make BLASLDFLAGS="${blasldflags}" ${flags} out/libscsdir.${dlext}
make BLASLDFLAGS="${blasldflags}" ${flags} out/libscsindir.${dlext}

mkdir -p ${libdir}
cp out/libscs*.${dlext} ${libdir}
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libscsindir", :libscsindir),
    LibraryProduct("libscsdir", :libscsdir)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("OpenBLAS_jll")
]

# Build the tarballs, and possibly a `build.jl` as well
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct BootDiagnosticsInstanceView <: SwaggerModel
    consoleScreenshotBlobUri::Any # spec type: Union{ Nothing, String } # spec name: consoleScreenshotBlobUri
    serialConsoleLogBlobUri::Any # spec type: Union{ Nothing, String } # spec name: serialConsoleLogBlobUri
    status::Any # spec type: Union{ Nothing, InstanceViewStatus } # spec name: status

    function BootDiagnosticsInstanceView(;consoleScreenshotBlobUri=nothing, serialConsoleLogBlobUri=nothing, status=nothing)
        o = new()
        validate_property(BootDiagnosticsInstanceView, Symbol("consoleScreenshotBlobUri"), consoleScreenshotBlobUri)
        setfield!(o, Symbol("consoleScreenshotBlobUri"), consoleScreenshotBlobUri)
        validate_property(BootDiagnosticsInstanceView, Symbol("serialConsoleLogBlobUri"), serialConsoleLogBlobUri)
        setfield!(o, Symbol("serialConsoleLogBlobUri"), serialConsoleLogBlobUri)
        validate_property(BootDiagnosticsInstanceView, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type BootDiagnosticsInstanceView

const _property_map_BootDiagnosticsInstanceView = Dict{Symbol,Symbol}(Symbol("consoleScreenshotBlobUri")=>Symbol("consoleScreenshotBlobUri"), Symbol("serialConsoleLogBlobUri")=>Symbol("serialConsoleLogBlobUri"), Symbol("status")=>Symbol("status"))
const _property_types_BootDiagnosticsInstanceView = Dict{Symbol,String}(Symbol("consoleScreenshotBlobUri")=>"String", Symbol("serialConsoleLogBlobUri")=>"String", Symbol("status")=>"InstanceViewStatus")
Base.propertynames(::Type{ BootDiagnosticsInstanceView }) = collect(keys(_property_map_BootDiagnosticsInstanceView))
Swagger.property_type(::Type{ BootDiagnosticsInstanceView }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_BootDiagnosticsInstanceView[name]))}
Swagger.field_name(::Type{ BootDiagnosticsInstanceView }, property_name::Symbol) =  _property_map_BootDiagnosticsInstanceView[property_name]

function check_required(o::BootDiagnosticsInstanceView)
    true
end

function validate_property(::Type{ BootDiagnosticsInstanceView }, name::Symbol, val)
end


text:
# Temporarily add this repo's Maple to the module path.
# Do not use this if you installed Maple using Pkg.clone.
push!(LOAD_PATH, "$(pwd())/src")

using Maple

# See top dictionary in `entity_map.jl` for currently implemented shorthands
lvl_1_fg = """
333333333333333333333333333333333333333333333333
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3                                              3
3333333                                        3
3vvvvvv    Q                         F         3
3                                              3
3        O                                     3
3>                                             3
3>    3333>  R                                 3
3     3333>        DDDDDD                      3
3                  DDDDDD                      3
3P ^^^^^^                                      3
333333333333333333333333333333333333333333333333
"""

fgTiles = FgTiles(lvl_1_fg)

map = Map(
    "EntityMapExample",
    Room[
        Room(
            name = "lvl_1",

            fgTiles = fgTiles,

            position = (0, 0),
            size = size(fgTiles),

            entities = entityMap(lvl_1_fg)
        )
    ]
)

encodeMap(map, "Testing-EntityMap.bin")

max_stars_repo_name: CelestialCartographers/Maple

max_issues_repo_name: CelestialCartographers/Maple

max_forks_repo_name: CelestialCartographers/Maple


text:
using StatsBase, PythonOT, Distances

export EarthMoversDistance, EMD
export check_trans_plan, get_info

# Optimal Transport (OT) Distances
# --------------------------------

abstract type LengthDistance <: Metric end

struct AbsoluteDiff <: LengthDistance end
function (d::AbsoluteDiff)(N::Int, M::Int)
    return abs(N-M)
end
struct SquaredDiff <: LengthDistance end
function (d::SquaredDiff)(N::Int, M::Int)
    return (N-M)^2
end

struct EarthMoversDistance{T<:Metric} <: Metric
    ground_dist::T
end

const EMD = EarthMoversDistance

function check_trans_plan(d::EMD, X::Vector{T}, Y::Vector{T}) where {T}

    a = proportionmap(X)
    b = proportionmap(Y)

    C = zeros(Float64, length(a), length(b))

    for (j, val_b) in enumerate(keys(b))
        for (i, val_a) in enumerate(keys(a))
            C[i, j] = d.ground_dist(val_a, val_b)
        end
    end

    TransPlan::Array{Float64,2} = PythonOT.emd(
        collect(values(a)),
        collect(values(b)), C
    )

    if sum(TransPlan) ≈ 1.0
        println("Transplan looks valid!")
    else
        error("Transplan looks spurious!")
    end
    return C
end

function (d::EMD)(X::Vector{T}, Y::Vector{T}) where {T}

    a = proportionmap(X)
    b = proportionmap(Y)

    C = zeros(Float64, length(a), length(b))

    for (j, val_b) in enumerate(keys(b))
        for (i, val_a) in enumerate(keys(a))
            C[i, j] = d.ground_dist(val_a, val_b)
        end
    end
    return PythonOT.emd2(
        collect(values(a)),
        collect(values(b)), C
    )
end

function get_info(
    d::EMD,
    X::Vector{T}, Y::Vector{T}
    ) where {T}

    a = proportionmap(X)
    b = proportionmap(Y)

    C = zeros(Float64, length(a), length(b))

    for (j, val_b) in enumerate(keys(b))
        for (i, val_a) in enumerate(keys(a))
            C[i, j] = d.ground_dist(val_a, val_b)
        end
    end

    TransPlan = PythonOT.emd(
        collect(values(a)),
        collect(values(b)), C
    )

    return collect(keys(a)), collect(keys(b)), TransPlan
end

# EMD composed with a distance of the number of interactions
struct sEMD{T<:Metric, G<:LengthDistance} <: Metric
    ground_dist::T
    length_dist::G
    τ::Real # Relative weighting term (a proportion weighting EMD vs length distance, high τ => high EMD weighting)
end

function (d::sEMD)(S1::Vector{T}, S2::Vector{T}) where {T}

    d₁ = EMD(d.ground_dist)(S1, S2)
    d₂ = d.length_dist(length(S1), length(S2))
    # @show d₁, d₂

    return d₁ + d.τ * d₂

end

struct sEMD2{T<:Metric, G<:LengthDistance} <: Metric
    ground_dist::T
    length_dist::G
    τ::Real
end

function (d::sEMD2)(S1::Vector{T}, S2::Vector{T}) where {T}

    d₁ = EMD(d.ground_dist)(S1, S2)
    d₂ = d.length_dist(sum(length.(S1)), sum(length.(S2)))
    # @show d₁, d₂

    return d.τ * d₁ + (1-d.τ) * d₂

end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
VERSION

# %%
indnearest_alloc(M, x) = argmin(@. abs(M - x)) # allocates!

# See https://github.com/JuliaLang/julia/blob/master/base/reduce.jl#L862
valindargmin(f, X) = mapfoldl(((i, x),) -> (f(x), i, x),
    ((v, i, x), (w, j, y)) -> v ≥ w ? (w, j, y) : (v, i, x), pairs(X))
valindargnearest(M, x) = valindargmin(m -> abs(m - x), M)
valindnearest(M, x) = ((val, ind, arg) = valindargnearest(M, x); (val, ind))
indargnearest(M, x) = ((val, ind, arg) = valindargnearest(M, x); (ind, arg))
indnearest(M, x) = ((val, ind, arg) = valindargnearest(M, x); ind)

# %%
f(x, y) = sin(x) * cos(x) + 0.01x * y
M = f.(2.5:0.2:3.7, (1:0.2:2)')

# %%
indnearest_alloc(M, 0.2)

# %%
Tuple(indnearest_alloc(M, 0.2))

# %%
valindargnearest(M, 0.2)

# %%
valindnearest(M, 0.2)

# %%
indargnearest(M, 0.2)

# %%
indnearest(M, 0.2)

# %%
Tuple(indnearest(M, 0.2))

# %%
using BenchmarkTools
@btime indnearest_alloc($M, 0.2)
@btime valindargnearest($M, 0.2)
@btime valindnearest($M, 0.2)
@btime indargnearest($M, 0.2)
@btime indnearest($M, 0.2)
@btime Tuple(indnearest($M, 0.2));

# %%

max_stars_repo_path: 0013/find nearest.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/find nearest.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/find nearest.jl

max_forks_repo_name: genkuroki/public


text:
# # Sparse matrices
#
# Source to source automatic differentiation is useful in differentiating sparse matrices. It is a well-known problem that sparse matrix operations can not benefit directly from generic backward rules for dense matrices because general rules do not keep the sparse structure.
# In the following, we will show that reversible AD can differentiate the Frobenius dot product between two sparse matrices with the state-of-the-art performance. Here, the Frobenius dot product is defined as \texttt{trace(A'B)}.
# Its native Julia (irreversible) implementation is `SparseArrays.dot`.
#
# The following is a reversible counterpart

using NiLang, NiLang.AD
using SparseArrays

@i function idot(r::T, A::SparseMatrixCSC{T},B::SparseMatrixCSC{T}) where {T}
    m ← size(A, 1)
    n ← size(A, 2)
    @invcheckoff branch_keeper ← zeros(Bool, 2*m)
    @safe size(B) == (m,n) || throw(DimensionMismatch("matrices must have the same dimensions"))
    @invcheckoff @inbounds for j = 1:n
        ia1 ← A.colptr[j]
        ib1 ← B.colptr[j]
        ia2 ← A.colptr[j+1]
        ib2 ← B.colptr[j+1]
        ia ← ia1
        ib ← ib1
        @inbounds for i=1:ia2-ia1+ib2-ib1-1
            ra ← A.rowval[ia]
            rb ← B.rowval[ib]
            if (ra == rb, ~)
                r += A.nzval[ia]' * B.nzval[ib]
            end
            ## b move -> true, a move -> false
            branch_keeper[i] ⊻= @const ia == ia2-1 || ra > rb
            ra → A.rowval[ia]
            rb → B.rowval[ib]
            if (branch_keeper[i], ~)
                INC(ib)
            else
                INC(ia)
            end
        end
        ~@inbounds for i=1:ia2-ia1+ib2-ib1-1
            ## b move -> true, a move -> false
            branch_keeper[i] ⊻= @const ia == ia2-1 || A.rowval[ia] > B.rowval[ib]
            if (branch_keeper[i], ~)
                INC(ib)
            else
                INC(ia)
            end
        end
    end
    @invcheckoff branch_keeper → zeros(Bool, 2*m)
end

# Here, the key point is using a \texttt{branch\_keeper} vector to cache branch decisions.

# The time used for a native implementation is

using BenchmarkTools
a = sprand(1000, 1000, 0.01);
b = sprand(1000, 1000, 0.01);
@benchmark SparseArrays.dot($a, $b)

# To compute the gradients, we wrap each matrix element with `GVar`, and send them to the reversible backward pass

out! = SparseArrays.dot(a, b)
@benchmark (~idot)($(GVar(out!, 1.0)),
        $(GVar.(a)), $(GVar.(b)))

# The time used for computing backward pass is approximately 1.6 times Julia's native forward pass.
# Here, we have turned off the reversibility check off to achieve better performance.
# By writing sparse matrix multiplication and other sparse matrix operations reversibly,
# we will have a differentiable sparse matrix library with proper performance.

# See my another blog post for [reversible sparse matrix multiplication](https://nextjournal.com/giggle/how-to-write-a-program-differentiably).


text:
# Autogenerated wrapper script for yaml_cpp_jll for x86_64-apple-darwin
export libyaml_cpp

JLLWrappers.@generate_wrapper_header("yaml_cpp")
JLLWrappers.@declare_library_product(libyaml_cpp, "@rpath/libyaml-cpp.0.7.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libyaml_cpp,
        "lib/libyaml-cpp.0.7.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""

`module Regularisers`

## Old Documentation Copy-Pasted ...

`regularise_2b(B, r0, r1; creg = 1e-2, Nquad = 20)`

construct a regularising stabilising matrix that acts only on 2-body
terms.

* `B` : basis
* `r0, r1` : upper and lower bound over which to integrate
* `creg` : multiplier
* `Nquad` : number of quadrature / sample points

```
   I = ∑_r h | ∑_j c_j ϕ_j''(r) |^2
     = ∑_{i,j} c_i c_j  ∑_r h ϕ_i''(r) ϕ_j''(r)
     = c' * M * c
   M_ij = ∑_r h ϕ_i''(r) ϕ_j''(r) = h * Φ'' * Φ''
   Φ''_ri = ϕ_i''(r)
```

`regularise(N, B, r0, r1; kwargs...)`

* `N::Integer` : body order
* `B::Vector` : basis
"""
module Regularisers

# TODO: move all this into NBodyIPs!

using StaticArrays
using JuLIP: AbstractCalculator
using JuLIP.Potentials: evaluate_d
using NBodyIPs: bodyorder, transform, evaluate_many_ricoords!, descriptor,
                  inv_transform, split_basis, BondLengthDesc, IdTransform,
                  BondAngleDesc
using NBodyIPs.Polys: NBPoly
using NBodyIPs.EnvIPs: EnvIP
using LinearAlgebra: I, Diagonal
using NBodyIPs.Sobol: filtered_sobol, filtered_cart_sobol, bl_is_simplex,
                      ba_is_simplex, Cart2BA

import Base: Matrix, Dict

export BLRegulariser, BLReg, BARegulariser, BAReg,
       Regulariser, EnvBLReg, EnvBLRegulariser,
       EnvBAReg, EnvBARegulariser

macro def(name, definition)
    return quote
        macro $(esc(name))()
            esc($(Expr(:quote, definition)))
        end
    end
end

abstract type NBodyRegulariser{N} end

abstract type EnvReg{N} <: NBodyRegulariser{N} end

@def nbregfields begin
   N::Int
   npoints::Int
   creg::T
   r0::T
   r1::T
   transform
   sequence::Symbol
   freg::Function
   valN::Val{N}
end

struct BLRegulariser{N, T} <: NBodyRegulariser{N}
   @nbregfields
   species::Vector{Int}
end

struct EnvBLRegulariser{N, T} <: EnvReg{N}
   @nbregfields
   envdeg::Int
   species::Vector{Int}
end

struct EnvBARegulariser{N, T} <: EnvReg{N}
   @nbregfields
   envdeg::Int
   species::Vector{Int}
end

struct BARegulariser{N, T} <: NBodyRegulariser{N}
   @nbregfields
   species::Vector{Int}
end

struct EnergyRegulariser{N, T} <: NBodyRegulariser{N}
   @nbregfields
   species::Vector{Int}
end

nbodyreg2dict(reg) = Dict(
   "type" => string(typeof(reg)),
   "N" => reg.N,
   "npoints" => reg.npoints,
   "creg" => reg.creg,
   "r0" => reg.r0,
   "r1" => reg.r1,
   "sequence" => string(reg.sequence),
   "transform" => Dict(reg.transform) )

Dict(reg::NBodyRegulariser) = nbodyreg2dict(reg)

function Dict(reg::EnvReg)
   D = nbodyreg2dict(reg)
   D["envdeg"] = reg.envdeg
end

const BLReg = BLRegulariser
const EnvBLReg = EnvBLRegulariser
const EnvBAReg = EnvBARegulariser
const BAReg = BARegulariser

BLRegulariser(N, r0, r1;
             creg = 1.0,
             npoints = Nquad(Val(N)),
             sequence = :sobol,
             transform = IdTransform(),
             freg = laplace_regulariser,
             species = []) =
   BLRegulariser(N, npoints, creg, r0, r1,
                 transform, sequence, freg, Val(N), species)

BARegulariser(N, r0, r1;
             npoints = Nquad(Val(N)),
             creg = 1.0,
             transform = IdTransform(),
             sequence = :sobol,
             freg = laplace_regulariser,
             species = []) =
   BARegulariser(N, npoints, creg, r0, r1,
                 transform, sequence, freg, Val(N), species)

EnvBLRegulariser(N, envdeg, r0, r1;
             creg = 1.0,
             npoints = Nquad(Val(N)),
             sequence = :sobol,
             transform = IdTransform(),
             freg = laplace_regulariser,
             species = [] ) =
   EnvBLRegulariser(N, npoints, creg, r0, r1,
                    transform, sequence, freg, Val(N), envdeg, species)

EnvBARegulariser(N, envdeg, r0, r1;
             creg = 1.0,
             npoints = Nquad(Val(N)),
             sequence = :sobol,
             transform = IdTransform(),
             freg = laplace_regulariser,
             species = []) =
   EnvBARegulariser(N, npoints, creg, r0, r1,
                    transform, sequence, freg, Val(N), envdeg, species)

# ===========================================================
#      originals ...

_bainvt(inv_t, x::StaticVector{1}) =
      inv_t.(x), SVector{0, Float64}()
_bainvt(inv_t, x::StaticVector{3}) =
      SVector(inv_t(x[1]), inv_t(x[2])), SVector(x[3])
_bainvt(inv_t, x::StaticVector{6}) =
      SVector(inv_t(x[1]), inv_t(x[2]), inv_t(x[3])),  SVector(x[4], x[5], x[6])
_bainvt(inv_t, x::StaticVector{10}) =
      SVector( inv_t(x[1]), inv_t(x[2]), inv_t(x[3]), inv_t(x[4]) ),
      SVector(x[5], x[6], x[7], x[8], x[9], x[10])

find_sub_basis(reg::NBodyRegulariser{N}, B) where {N} =
   findall(b -> ( (bodyorder(b) .== N) &&
          ((reg.species == []) || (species(b) == reg.species))
          ), B)

# find_sub_basis(reg::EnvReg{N}, B) where {N} =
#    findall( b -> ( (bodyorder(b) == N) && (b.t == reg.envdeg) ), B )

find_sub_basis(reg::EnvReg{N}, B) where {N} =
   findall( b -> ( (bodyorder(b) == N) &&
                   (b.t == reg.envdeg) &&
                   (  (reg.species == []) || (species(b) == reg.species) )
           ), B )

species(b::NBPoly) = []

#
# this converts the Regulariser type information to a matrix that can be
# attached to the LSQ problem.
#
function Matrix(reg::NBodyRegulariser{N}, B::Vector{<: AbstractCalculator};
                verbose=false
                ) where {N}

   # TODO: let the regulariser decide which basis functions it can
   #       be applied to, or let the user adjust it!
   Ib = find_sub_basis(reg, B)
   if isempty(Ib)
      verbose && @warn("""Trying to construct a $N-body regulariser, but no basis
                         function with bodyorder $N exists.""")
      Ψreg = zeros(0, length(B))
      Yreg = zeros(0, 1)
      return Ψreg, Yreg
   end
   # scalar inverse transform
   inv_t = x -> inv_transform(reg.transform, x)

   # filter
   if reg isa Union{BLRegulariser, EnvBLRegulariser}
      # vectorial inverse transform
      inv_tv = x -> inv_t.(x)
      filter = x -> bl_is_simplex( inv_tv(x) )
      x0 = transform(reg.transform, reg.r0) * SVector(ones((N*(N-1))÷2)...)
      x1 = transform(reg.transform, reg.r1) * SVector(ones((N*(N-1))÷2)...)
   elseif reg isa Union{BARegulariser, EnvBARegulariser}
      inv_tv = x -> _bainvt(inv_t, x)
      filter = x -> ba_is_simplex( inv_tv(x)... )
      o1 = SVector(ones(N-1)...)
      o2 = N > 2 ? SVector(ones( ((N-1)*(N-2))÷2 )...) : SVector{0,Float64}()
      x0 = vcat( transform(reg.transform, reg.r0) * o1, -o2 )
      x1 = vcat( transform(reg.transform, reg.r1) * o1, o2 )
   else
      @show reg
      @error("Unknown type of reg: `typeof(reg) == $(typeof(reg))`")
   end

   # choose the sub-basis to evaluate, normally that's just B[Ib], but
   #
   if reg isa Union{BLRegulariser, BARegulariser}
      # we do this in a weird way since the array needs to infer the type
      # but B[Ib] doesn't for some reason.
      subB = [b for b in B[Ib]]
   elseif reg isa Union{EnvBLRegulariser, EnvBARegulariser}
      subB = [b.Vr for b in B[Ib]]
   else
      @error("Unknown type of reg")
   end

   if N == 2
      # uniform grid in 2d...
      X = SVector.(range(x0, x1, length = reg.npoints))
   else
      if reg.sequence == :sobol
         # construct a low discrepancy sequence
         X = filtered_sobol(x0, x1, filter; npoints=reg.npoints, nfailed=100*reg.npoints)
      elseif reg.sequence == :cartsobol
         # construct a low discrepancy sequence
         @assert reg isa BARegulariser
         converter = Cart2BA(N)
         filter = x -> all( (x0[i] < x[1][i] < x1[i]) for i = 1:N-1 )
         X = filtered_cart_sobol(maximum(x1), length(x1), converter, filter,
                                 npoints = reg.npoints,
                                 nfailed = 1000*reg.npoints)
      elseif reg.sequence == :cart
         error("TODO: implement `sequence == :cart`")
      elseif reg.sequence isa Vector
         # this is intended to put the regulariser on configs in the database
         error("TODO: implement `sequence isa Vector`")
      else
         error("unknown argument `sequence = $sequence`")
      end
   end

   # loop through sobol points and collect the laplacians at each point.
   Ψreg = reg.creg/length(X) * assemble_reg_matrix(X, subB, length(B),
                                                   Ib, inv_tv, reg.freg)
   Yreg = zeros(size(Ψreg, 1))
   return Ψreg, Yreg
end

_envdeg_(b::NBPoly) = 0
_envdeg_(b::EnvIP) = b.t

regeval_d(b::NBPoly, args...) = evaluate_d(b, args...)
regeval_d(b::EnvIP, args...) = evaluate_d(b.Vr, args...)

_Vr(b::NBPoly) = b
_Vr(b::EnvIP) = b.Vr

# X = Sobol points
# B = basis[Ib]
# nB = length(basis)
# inv_tv = inverse space transform
# freg = regularisation function (e.g. laplace)
#
function assemble_reg_matrix(X, B, nB, Ib, inv_tv, freg)
   @assert length(B) == length(Ib)
   # split the basis into "nice" parts
   Bord, Iord = split_basis(B; splitfun = b -> typeof(b))
   # allocate regularisation matrices for each of these
   Ψ = zeros(length(X), nB)
   # go through the basis subsets
   for n = 1:length(Bord)
      # indices of basis subset Bord[n] in basis
      Jb = Ib[Iord[n]]
      # temprary storage for the regularisation
      temp = zeros(length(Jb))
      # loop through points x ∈ X at which to apply the regularisation
      # function freg
      for (ix, x) in enumerate(X)
         Ψ[ix, Jb] = freg(x, Bord[n], temp, inv_tv)
      end
   end
   return Ψ
end

function laplace_regulariser(x::SVector{DIM,T}, B::Vector{<: AbstractCalculator},
                             temp::Vector{T}, inv_tv) where {DIM, T}
   if !isconcretetype(eltype(B))
      @warn("laplace_regulariser: `TB` is not a leaf type")
   end
   h = 1e-2
   rθ = inv_tv(x)
   evaluate_many_ricoords!(temp, B, rθ)
   L = 2 * DIM * copy(temp)

   EE = SMatrix{DIM,DIM}(1.0I)
   for j = 1:DIM
      rθp = inv_tv(x + h * EE[:,j])
      evaluate_many_ricoords!(temp, B, rθp)
      L -= temp
      rθm = inv_tv(x - h * EE[:,j])
      evaluate_many_ricoords!(temp, B, rθm)
      L -= temp
   end

   return L/h^2
end

# laplace_regulariser(::Tuple{StaticArrays.SArray{Tuple{5},Float64,1,5},StaticArrays.SArray{Tuple{10},Float64,1,10}},
#                     ::Array{NBPoly{5,11,Float64,BondAngleDesc{PolyTransform{Int64,Float64},PolyCut2sA{Float64}}},1},
#                     ::Array{Float64,1},
#                     ::getfield(NBodyIPs.Regularisers, Symbol("##11#21"))
#                  )
# ==================================================================
#               ENERGY REGULARISER
# ==================================================================

# TODO: revisit this one!

# EnergyRegulariser(N, r0, r1;
#             npoints = Nquad(Val(N)),
#             creg = 0.1,
#             sequence = :sobol) =
#    EnergyRegulariser(N, npoints, creg, r0, r1, sequence, energy_regulariser, Val(N))

# function energy_regulariser(x::SVector{DIM,T}, B::Vector{<: AbstractCalculator},
#                             temp::Vector{T}, inv_tv) where {DIM, T}
#    evaluate_many_ricoords!(temp, B, r)
#    return - temp
# end

# ==================================================================
#               DEFAULTS
# ==================================================================

Nquad(::Val{2}) = 100
Nquad(::Val{3}) = 1000
Nquad(::Val{4}) = 10_000

Regulariser(V::NBPoly{N, M, T, TD}, r0, r1; kwargs...
           ) where {N, M, T} where {TD <: BondLengthDesc} =
   BLReg(N, r0, r1; creg=1.0, npoints = Nquad(Val(N)),
                    transform = transform(V), kwargs... )

Regulariser(V::NBPoly{N, M, T, TD}, r0, r1; kwargs...
           ) where {N, M, T} where {TD <: BondAngleDesc} =
   BAReg(N, r0, r1; creg=1.0, npoints = Nquad(Val(N)),
                    transform = transform(V), kwargs... )

# ==================================================================

"""
`struct L2Regulariser`: standard Tikhonov regulariser

## Constructors

```
L2Regulariser(α::Float64) # => α² Σᵢ xᵢ²
L2Regulariser(a::Vector)  # => Σᵢ aᵢ² xᵢ²
L2Regulariser(db, A::Dict)  # see below
```
The first constructor is the naive Tikhonov constructor with a constant
regularisation parameter.  The second constructor gives maximal flexibility
allowing prescribing an individual parameter to any basis function.
In the last constructor A[n] is the Tikhonov parameter for basis
functions withbody-order n. The lsq system stored in `db` is scanned
for indices of these basis functions and a vector `a` is  then formed
to call the second constructor.
"""
struct L2Regulariser
   a::Vector{Float64}
   α::Float64
end

L2Regulariser(α::Float64) = L2Regulariser(Float64[], α)

L2Regulariser(a::Vector) = L2Regulariser(a, 0.0)

function L2Regulariser(db, A::Dict{Int, <: Real})
   a = zeros(length(db.basis))
   for (n, α) in A
      In = findall(b -> bodyorder(b) == n, db.basis)
      a[In] .= α
   end
   return L2Regulariser(a)
end

function Matrix(reg::L2Regulariser, B::Vector{<: AbstractCalculator}; kwargs...)
   if isempty(reg.a)
      return Matrix( a*I, (length(B), length(B)) ), zeros(length(B))
   else
      return collect(Diagonal(reg.a)), zeros(length(B))
   end
end

Dict(reg::L2Regulariser) = Dict("__id__" => "L2Regulariser",
                                "a" => reg.a, "α" => reg.α)

end

max_stars_repo_name: cortner/ManyBodyIPs

max_issues_repo_name: cortner/ManyBodyIPs

max_forks_repo_name: cortner/ManyBodyIPs


text:
include("rendering.jl")

mutable struct Lightray
    x::Float64
    y::Float64
    angle::Float64
end

mutable struct Tree
    x::Float64
    y::Float64
    growth::Float64
end

function drawlightray(r, lightray)
    drawsetcolor(r, 255, 255, 255)
    x2 = lightray.x + cos(lightray.angle) * 20
    y2 = lightray.y - sin(lightray.angle) * 20
    drawline(r, lightray.x, lightray.y, x2, y2)
end

function drawtree(r, x, y, size, depth, angle, setofrays)
    if (depth > 0)
        drawsetcolor(r, 255, 255, 0)
        angleinc = 0.6
        langle = angle + angleinc
        rangle = angle - angleinc
        eplx = x - cos(langle)*size
        eply = y - sin(langle)*size
        eprx = x - cos(rangle)*size
        epry = y - sin(rangle)*size
        drawline(r,x,y,eplx,eply)
        drawline(r,x,y,eprx,epry)
        accum = 0
        accum = accum + drawtree(r,eplx,eply,size/1.6,depth-1,langle,setofrays)
        accum = accum + drawtree(r,eprx,epry,size/1.6,depth-1,rangle,setofrays)
        return accum
    else
        drawsetcolor(r, 0, 255, 0)
        drawline(r,x-10,y-10,x+10,y-10)
        drawline(r,x-10,y+10,x+10,y+10)
        drawline(r,x-10,y-10,x-10,y+10)
        drawline(r,x+10,y-10,x+10,y+10)
        accum = 0
        todelete = Set()
        for ray in setofrays
            if (ray.x > x-10 && ray.x < x + 10 && ray.y > y-10 && ray.y < y+10)
                accum = accum + 1
                push!(todelete,ray)
            end
        end
        for ray in todelete
            delete!(setofrays,ray)
        end
        return accum
    end
end

function lerp(a,b,progress)
    initialresult = (b-a)*progress
    return a+initialresult
end

function createrays!(ox,oy,set,angle,size)
    for i in 1:2
        epointx = ox + cos(angle) * size
        epointy = oy - sin(angle) * size
        lpointx = epointx + cos(angle+pi/2) * size
        lpointy = epointy - sin(angle+pi/2) * size
        rpointx = epointx + cos(angle-pi/2) * size
        rpointy = epointy - sin(angle-pi/2) * size
        progress = rand(Float32)
        xfinal = lerp(lpointx, rpointx, progress)
        yfinal = lerp(lpointy, rpointy, progress)
        rayangle = angle + pi
        push!(set,Lightray(xfinal,yfinal,rayangle))
    end
end

function updaterays!(setofrays,treebaseline,lox,loy,lsize,progress,r)
    createrays!(400,400,setofrays,progress,600)
    for thisray in setofrays
        drawlightray(r,thisray)
        thisray.x = thisray.x + cos(thisray.angle) * 5
        thisray.y = thisray.y - sin(thisray.angle) * 5
        if thisray.y > treebaseline
            delete!(setofrays,thisray)
        end
        if thisray.x > lox + lsize + 100
            delete!(setofrays,thisray)
        end
        if thisray.x < lox - lsize - 100
            delete!(setofrays,thisray)
        end
        if thisray.y < loy - lsize - 100
            delete!(setofrays,thisray)
        end
    end
end

function updatetree!(tree,r,setofrays)
    size = lerp(40,80,tree.growth)
    depth = trunc(Int,lerp(2,5,tree.growth))
    if depth > 6
        depth = 6
    end
    tree.growth = tree.growth + drawtree(r,tree.x,tree.y,size,depth,pi/2,setofrays)/(100*(tree.growth+0.1))
    tree.growth = tree.growth - 0.0005 * tree.growth
    if tree.growth < 0
        tree.growth = 0
    end
end

function treelightsim()
    r = drawinit()
    progress = 0
    setofrays = Set()
    treebaseline = 700
    tree1 = Tree(300, treebaseline, 0)
    tree2 = Tree(400, treebaseline, 0)
    lox = 400
    loy = 400
    lsize = 600
    function todraw()
        for i in 1:4
            progress += 0.001
            drawsetcolor(r, 0, 0, 0)
            drawclear(r)
            drawsetcolor(r, 255, 255, 255)
            updatetree!(tree1,r,setofrays)
            updatetree!(tree2,r,setofrays)
            updaterays!(setofrays,treebaseline,lox,loy,lsize,progress,r)
        end
    end
    drawloop(r, todraw)
end


text:
using MatrixCoalescent
using Test

@testset "Rational" begin
    n=3
    mc = MatCoal(Rational{Int}, n)
    @test mc.nLineages == 3
    @test dim(mc) == n-1
    @test mc.beta[1] == 1//1
    @test mc.beta[2] == 3//1

    @test mc.rbeta[1] == 1//1
    @test mc.rbeta[2] == 1//3

    @test mc.gmat[1,1] == 3//2
    @test mc.gmat[1,2] == -3//2
    @test mc.gmat[2,1] == 0//1
    @test mc.gmat[2,2] == 1//1

    @test mc.hmat[1,1] == -3//2
    @test mc.hmat[1,2] == 1//2
    @test mc.hmat[2,1] == 0//1
    @test mc.hmat[2,2] == -1//3
end

@testset "Float64" begin
    n=3
    mc = MatCoal(Float64, n)
    @test mc.nLineages == 3
    @test dim(mc) == n-1
    @test mc.beta[1] == 1.0
    @test mc.beta[2] == 3.0

    @test mc.rbeta[1] == 1.0
    @test mc.rbeta[2] == 1.0/3.0

    @test mc.gmat[1,1] == 1.5
    @test mc.gmat[1,2] == -3.0/2.0
    @test mc.gmat[2,1] == 0.0/1.0
    @test mc.gmat[2,2] == 1.0

    @test mc.hmat[1,1] == -1.5
    @test mc.hmat[1,2] == 0.5
    @test mc.hmat[2,1] == 0.0
    @test mc.hmat[2,2] == -1.0/3.0

    eig = Vector{Float64}(undef, n-1)
    v = 1.0
    eigenvals!(eig, v, mc)

    @test eig[1] == exp(-1)
    @test eig[2] == exp(-3)

    pr = Vector{Float64}(undef, n-1)
    project!(pr, eig, mc)

    @test pr[1] ≈ 1.5 * (exp(-1) - exp(-3))
    @test pr[2] ≈ exp(-3)

    elen = Vector{Float64}(undef, n-1)
    interval_lengths!(elen, eig, mc)

    elen[1] ≈ 1.0 - 1.5*exp(-1) + 0.5*exp(-3)
    elen[2] ≈ 1.0/3.0 - exp(-3)/3.0
end

"""
For a given number of lineages (nLin) and across a range of time values,
calculate the projected probability vector (using project!) and the
expected interval lengths (using interval_lengths!) using 64-bit
floats and also using 256-bit floats. Return the largest absolute
difference between the low-precision and the high-precision versions
of this calculation.
"""
function estimate_error(nLin::Int)
    mcl = MatCoal(Float64, nLin)  # low precision
    mch = MatCoal(BigFloat, nLin) # high precision

    el = Vector{Float64}(undef, nLin-1)
    eh = Vector{BigFloat}(undef, nLin-1)
    xl = Vector{Float64}(undef, nLin-1)
    xh = Vector{BigFloat}(undef, nLin-1)

    err = 0.0

    for v in 0.0 : 0.5 : 9.5
        eigenvals!(el, v, mcl)
        eigenvals!(eh, BigFloat(v), mch)
        project!(xl, el, mcl)
        project!(xh, eh, mch)
        for (l,h) in zip(xl, xh)
            err = max(err, abs(l - Float64(h)))
        end
        interval_lengths!(xl, el, mcl)
        interval_lengths!(xh, eh, mch)
        for (l,h) in zip(xl, xh)
            err = max(err, abs(l - Float64(h)))
        end
    end
    return err
end

@testset "Accuracy" begin
    @test estimate_error(2)  < eps(Float64)
    @test estimate_error(4)  < eps(Float64)
    @test estimate_error(8)  < 4e-15
    @test estimate_error(16) < 3.5e-13
    @test estimate_error(32) < 4e-8
    @test estimate_error(34) < 6e-8
    @test estimate_error(35) < 3e-7
    @test_throws OverflowError estimate_error(36)
end


text:
# this produces critical diagrams, however the rank values are input by hand,
# so thos needs to have some more tuning
using DrWatson
@quickactivate
using GenerativeAD

# save path
savepath = datadir("evaluation/cd_diagrams")
mkpath(savepath)

# number of datasets
nd = 40

# tabular
model_names = ["aae", "avae", "gano", "vae", "wae", "abod", "hbos", "if", "knn",
"loda", "lof", "orbf", "osvm", "pidf", "maf", "rnvp", "sptn", "fmgn", "gan",
"mgal", "dsvd", "vaek", "vaeo"]
model_names_gray = ["aae", "avae", "gano", "vae", "wae", "\\textcolor{gray}{abod}",
"\\textcolor{gray}{hbos}", "\\textcolor{gray}{if}", "\\textcolor{gray}{knn}",
"\\textcolor{gray}{loda}", "\\textcolor{gray}{lof}", "\\textcolor{gray}{orbf}",
"\\textcolor{gray}{osvm}", "\\textcolor{gray}{pidf}", "maf", "rnvp", "sptn", "fmgn",
"gan", "mgal", "dsvd", "vaek", "vaeo"]
model_ranks = [8.8, 13.4, 10.2, 8.0, 6.3, 11.3, 14.4, 14.2, 8.0, 15.6, 11.3, 8.1, 2.8,
	14.1, 8.8, 8.7, 10.3, 11.3, 11.8, 21.4, 16.0, 10.9, 6.7]
plot5 = GenerativeAD.Evaluation.ranks2tikzcd(model_ranks, model_names_gray,
	GenerativeAD.Evaluation.nemenyi_cd(length(model_ranks), nd, 0.05);
	scale=0.45, vpad=0.6, levelwidth=0.03, vdist=0.15)
GenerativeAD.Evaluation.string2file(joinpath(savepath, "tabular_auc_auc_cdd5.tex"), plot5)

plot10 = GenerativeAD.Evaluation.ranks2tikzcd(model_ranks, model_names_gray,
	GenerativeAD.Evaluation.nemenyi_cd(length(model_ranks), nd, 0.1);
	scale=0.45, vpad=0.6, levelwidth=0.03, vdist=0.15)
GenerativeAD.Evaluation.string2file(joinpath(savepath, "tabular_auc_auc_cdd10.tex"), plot10)

#
model_ranks_tpr5 = [10.5, 12.5, 9.9, 7.6, 8.2, 12.0, 14.7, 15.9, 9.8, 17.9, 10.9, 10.2,
	3.0, 13.9, 10.0, 9.8, 11.8, 9.6, 10.0, 18.9, 15.9, 11.2, 6.4]
plot5tpr = GenerativeAD.Evaluation.ranks2tikzcd(model_ranks_tpr5, model_names_gray,
	GenerativeAD.Evaluation.nemenyi_cd(length(model_ranks), nd, 0.05);
	scale=0.55, vpad=0.5, levelwidth=0.03, vdist=0.13)
GenerativeAD.Evaluation.string2file(joinpath(savepath, "tabular_tpr5_tpr5_cdd5.tex"),
	plot5tpr)

plot10tpr = GenerativeAD.Evaluation.ranks2tikzcd(model_ranks_tpr5, model_names_gray,
	GenerativeAD.Evaluation.nemenyi_cd(length(model_ranks), nd, 0.1);
	scale=0.55, vpad=0.5, levelwidth=0.03, vdist=0.13)
GenerativeAD.Evaluation.string2file(joinpath(savepath, "tabular_tpr5_tpr5_cdd10.tex"),
	plot10tpr)

# image data
nd = 40
img_model_names = ["aae", "gano", "skip", "vae", "wae", "\\textcolor{gray}{knn}",
"\\textcolor{gray}{osvm}", "fano", "fmgn", "dsvd", "vaek", "vaeo"]
img_model_ranks = [3.0, 6.2, 11.8, 6.0, 8.2, 8.0, 7.2, 5.8, 3.0, 6.2, 5.0, 2.5]
plot5img = GenerativeAD.Evaluation.ranks2tikzcd(img_model_ranks,img_model_names,
	GenerativeAD.Evaluation.nemenyi_cd(length(img_model_ranks), nd, 0.05))
GenerativeAD.Evaluation.string2file(joinpath(savepath, "images_auc_auc_cdd5.tex"),
	plot5img)

plot10img = GenerativeAD.Evaluation.ranks2tikzcd(img_model_ranks,img_model_names,
	GenerativeAD.Evaluation.nemenyi_cd(length(img_model_ranks), nd, 0.1))
GenerativeAD.Evaluation.string2file(joinpath(savepath, "images_auc_auc_cdd10.tex"),
	plot10img)

plot10img = GenerativeAD.Evaluation.ranks2tikzcd(img_model_ranks,img_model_names,
	GenerativeAD.Evaluation.nemenyi_cd(length(img_model_ranks), 4, 0.1))
GenerativeAD.Evaluation.string2file(
	joinpath(savepath, "images_auc_auc_cdd10_4_datasets.tex"), plot10img)

# img data per class
nd = 40
img_model_names = ["aae", "gano", "skip", "vae", "wae", "\\textcolor{gray}{knn}",
"\\textcolor{gray}{osvm}", "fano", "fmgn", "dsvd", "vaek", "vaeo"]
img_model_ranks_per_class = [5.0, 6.8, 10.6, 5.8, 6.2, 6.3, 7.4, 6.4,
	4.6, 6.5, 3.9, 3.2]
plot10img = GenerativeAD.Evaluation.ranks2tikzcd(img_model_ranks_per_class,img_model_names,
	GenerativeAD.Evaluation.nemenyi_cd(length(img_model_ranks_per_class), nd, 0.1))
GenerativeAD.Evaluation.string2file(joinpath(savepath, "images_auc_auc_cdd10_per_class.tex"),
	plot10img)


text:
using Test

@testset "Hyperslice solver" begin include("TestSolver.jl") end
@testset "Simplex/point intersection" begin include("TestSimplexPointIntersection.jl") end
@testset "Additional focus point sampling" begin include("TestSampling.jl") end
@testset "Utilities" begin include("TestUtil.jl") end
@testset "Top-level hyperslice function" begin include("TestHyperslice.jl") end

max_stars_repo_name: gabysbrain/multid-slicing

max_issues_repo_name: gabysbrain/multid-slicing


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct ManagementPoliciesApi <: SwaggerApi
    client::Swagger.Client
end

"""

Sets the managementpolicy to the specified storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: managementPolicyName::String (required)
Param: properties::ManagementPolicy (required)
Return: ManagementPolicy
"""
function _swaggerinternal_managementPoliciesCreateOrUpdate(_api::ManagementPoliciesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String, properties; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "managementPoliciesCreateOrUpdate", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "managementPoliciesCreateOrUpdate", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "managementPoliciesCreateOrUpdate", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "managementPoliciesCreateOrUpdate", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "managementPoliciesCreateOrUpdate", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "managementPoliciesCreateOrUpdate", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "PUT", ManagementPolicy, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/managementPolicies/{managementPolicyName}", ["azure_auth"], properties)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "managementPolicyName", managementPolicyName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function managementPoliciesCreateOrUpdate(_api::ManagementPoliciesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String, properties; _mediaType=nothing)
    _ctx = _swaggerinternal_managementPoliciesCreateOrUpdate(_api, resourceGroupName, accountName, api_version, subscriptionId, managementPolicyName, properties; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function managementPoliciesCreateOrUpdate(_api::ManagementPoliciesApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String, properties; _mediaType=nothing)
    _ctx = _swaggerinternal_managementPoliciesCreateOrUpdate(_api, resourceGroupName, accountName, api_version, subscriptionId, managementPolicyName, properties; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the managementpolicy associated with the specified storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: managementPolicyName::String (required)
Return: Nothing
"""
function _swaggerinternal_managementPoliciesDelete(_api::ManagementPoliciesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "managementPoliciesDelete", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "managementPoliciesDelete", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "managementPoliciesDelete", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "managementPoliciesDelete", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "managementPoliciesDelete", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "managementPoliciesDelete", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/managementPolicies/{managementPolicyName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "managementPolicyName", managementPolicyName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function managementPoliciesDelete(_api::ManagementPoliciesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_managementPoliciesDelete(_api, resourceGroupName, accountName, api_version, subscriptionId, managementPolicyName; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function managementPoliciesDelete(_api::ManagementPoliciesApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_managementPoliciesDelete(_api, resourceGroupName, accountName, api_version, subscriptionId, managementPolicyName; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the managementpolicy associated with the specified storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: managementPolicyName::String (required)
Return: ManagementPolicy
"""
function _swaggerinternal_managementPoliciesGet(_api::ManagementPoliciesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "managementPoliciesGet", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "managementPoliciesGet", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "managementPoliciesGet", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "managementPoliciesGet", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "managementPoliciesGet", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "managementPoliciesGet", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", ManagementPolicy, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/managementPolicies/{managementPolicyName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "managementPolicyName", managementPolicyName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function managementPoliciesGet(_api::ManagementPoliciesApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_managementPoliciesGet(_api, resourceGroupName, accountName, api_version, subscriptionId, managementPolicyName; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function managementPoliciesGet(_api::ManagementPoliciesApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, managementPolicyName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_managementPoliciesGet(_api, resourceGroupName, accountName, api_version, subscriptionId, managementPolicyName; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export managementPoliciesCreateOrUpdate, managementPoliciesDelete, managementPoliciesGet