text:
"""
$(TYPEDEF)

A triangular or symmetric matrix which stores data non-redundantly in a
contiguous linear array.
"""
struct TriMatrix{T, L<:TriLayout, A} <: AbstractMatrix{T}
	n::Int
	data::A
	diag::T

	"""
		TriMatrix(layout::TriLayout, n::Integer, data::AbstractVector{T}; diag=zero(T))

	Create an `n` by `n` `TriMatrix` using the existing 1D array `data` for storage.

	The length of `data` must match `n`, use [`TriMatrices.nelems`](@ref) to determine this value.
	"""
	function TriMatrix(layout::TriLayout, n::Integer, data::AbstractVector; diag=zero(eltype(data)))
		Base.require_one_based_indexing(data)
		length(data) == nelems(layout, n) || error("Matrix size does not match size of data array")
		return new{eltype(data), typeof(layout), typeof(data)}(n, data, diag)
	end
end

const TriUpperMatrix{T, D, A} = TriMatrix{T, TriUpper{D}, A}
const TriLowerMatrix{T, D, A} = TriMatrix{T, TriLower{D}, A}
const TriULMatrix{T, D, A} = TriMatrix{T, <:TriUL{D}, A}
const TriSymmetricMatrix{T, D, A} = TriMatrix{T, TriSymmetric{D}, A}

########################################
# Constructors
########################################

# Construct uninitialized
"""
	TriMatrix{T}(layout::TriLayout, undef, n::Integer; diag=zero(T))
	TriMatrix(layout::TriLayout, undef, n::Integer; diag=0.)

Create an uninitialized `n` by `n` `TriMatrix`.

`T` defaults to `Float64` if not specified.
"""
function TriMatrix{T}(layout::TriLayout, ::UndefInitializer, n::Integer; diag=zero(T)) where T
	data = Array{T}(undef, nelems(layout, n))
	return TriMatrix(layout, n, data, diag=convert(T, diag))
end

TriMatrix(layout::TriLayout, ::UndefInitializer, n::Integer; kw...) = TriMatrix{Float64}(layout, undef, n; kw...)

"""
	TriMatrix{T}(layout::TriLayout, m::AbstractMatrix; diag=zero(T))
	TriMatrix(layout::TriLayout, m::AbstractMatrix; [diag])

Create a `TriMatrix` by copying data from a square matrix `m`.

If unspecified `T` defaults to `eltype(m)`.
"""
function TriMatrix{T}(layout::TriLayout, m::AbstractMatrix;
                      diag=isempty(m) ? zero(T) : convert(T, m[1, 1])) where T
	n = LinearAlgebra.checksquare(m)

	tm = TriMatrix{T}(layout, undef, n, diag=diag)
	copyto!(tm, m)

	return tm
end

TriMatrix(layout::TriLayout, m::AbstractMatrix{T}; kw...) where T = TriMatrix{T}(layout, m; kw...)

"""
	fill(x, layout::TriLayout, n::Integer; diag=x)

Create an `n` by `n` [`TriMatrix`](@ref) with the stored data region filled with the value `x`.
"""
function Base.fill(x, layout::TriLayout, n::Integer; diag=x)
	mat = TriMatrix{typeof(x)}(layout, undef, n; diag=diag)
	fill!(mat.data, x)
	return mat
end

"""
	ones(T::Type=Float64, layout::TriLayout, n::Integer; diag=x)

Create an `n` by `n` [`TriMatrix`](@ref) with the stored data region filled with ones.
"""
Base.ones(T::Type, layout::TriLayout, n::Integer; diag=one(T)) = fill(one(T), layout, n, diag=diag)
Base.ones(layout::TriLayout, n::Integer; kw...) = ones(Float64, layout, n; kw...)

"""
	zeros(T::Type=Float64, layout::TriLayout, n::Integer; diag=x)

Create an `n` by `n` [`TriMatrix`](@ref) with the stored data region filled with zeros.
"""
Base.zeros(T::Type, layout::TriLayout, n::Integer; diag=zero(T)) = fill(zero(T), layout, n, diag=diag)
Base.zeros(layout::TriLayout, n::Integer; kw...) = zeros(Float64, layout, n; kw...)

# Construct similar
Base.similar(m::TriMatrix, new_eltype::Type, dims::Dims) = similar(m.data, new_eltype, dims)

########################################
# General methods
########################################

"""
	TriMatrices.TriLayout(m::TriMatrix)

Get the `TriLayout` of the given `TriMatrix` instance.
"""
TriLayout(m::TriMatrix) = TriLayout(typeof(m))
TriLayout(::Type{<:TriMatrix{T, L}}) where {T, L} = L()

TriMatrices.Indexing.tri_indices(m::TriMatrix) = tri_indices(TriLayout(m), m.n)

Base.size(m::TriMatrix) = (m.n, m.n)
Base.IndexStyle(::Type{<:TriMatrix}) = IndexCartesian()

"""
	parent(m::TriMatrix)

Get the underlying data array of the [`TriMatrix`](@ref) `m`.
"""
Base.parent(m::TriMatrix) = m.data

# Copy data from general matrix into TriMatrix
function Base.copyto!(dest::TriMatrix, src::AbstractMatrix)
	n = dest.n
	size(src) == (n, n) || error("Matrix sizes do not match")

	for (i, (r, c)) in enumerate(tri_indices(dest))
		@inbounds dest.data[i] = src[r, c]
	end

	return dest
end

########################################
# Indexing
########################################

"""
$(TYPEDSIGNATURES)

Get the value of an element within the stored data region of a
[`TriMatrix`](@ref). Should be faster than `Base.getindex` (even with
`@inbounds` applied) but results are undefined if `check_tri_index(mat, i, j)`
does not pass.
"""
function getindex_tri_unsafe(mat::TriMatrix, i::Integer, j::Integer)
	ii = car2lin_unchecked(TriLayout(mat), i, j)
	return @inbounds mat.data[ii]
end
getindex_tri_unsafe(mat::TriMatrix, idx::CartesianIndex{2}) = getindex_tri_unsafe(mat, idx[1], idx[2])

@Base.propagate_inbounds function Base.getindex(mat::TriMatrix{T,L},
                                                i::Integer, j::Integer,
                                                ) where {T, L}

	@boundscheck checkbounds(mat, i, j)

	!hasdiag(L) && i == j && return mat.diag  # Diagonal fill value
	L <: TriLower && i < j && return zero(T)  # Lower triangular above the diagonal
	L <: TriUpper && i > j && return zero(T)  # Upper triangular below the diagonal

	return getindex_tri_unsafe(mat, i, j)
end

"""
$(TYPEDSIGNATURES)

Set the value of an element within the stored data region of a
[`TriMatrix`](@ref). Should be faster than `Base.setindex!` (even with
`@inbounds` applied) but results are undefined if `check_tri_index(mat, i, j)`
does not pass.
"""
function setindex_tri_unsafe!(mat::TriMatrix, v, i::Integer, j::Integer)
	ii = car2lin_unchecked(TriLayout(mat), i, j)
	return @inbounds mat.data[ii] = v
end
setindex_tri_unsafe!(mat::TriMatrix, v, idx::CartesianIndex{2}) = setindex_tri_unsafe!(mat, v, idx[1], idx[2])

@Base.propagate_inbounds function Base.setindex!(mat::TriMatrix{T,L}, v,
                                                 i::Integer, j::Integer,
                                                ) where {T, L}

	@boundscheck checkbounds(mat, i, j)

	if !hasdiag(L) && i == j
		v != mat.diag && error("Cannot change value on diagonal for TriMatrix with layout $L")
		return mat.diag
	end

	if L <: TriLower && i < j
		v != 0 && error("Cannot assign non-zero value to index above diagonal for TriMatrix with layout TriLower")
		return zero(T)
	end

	if L <: TriUpper && i > j
		v != 0 && error("Cannot assign non-zero value to index below diagonal for TriMatrix with layout TriUpper")
		return zero(T)
	end

	return setindex_tri_unsafe!(mat, v, i, j)
end

########################################
# Display
########################################

# Renders entries on zero side of diagonal as dots, like triangular types in LinearAlgebra
function Base.replace_in_print_matrix(::TriMatrix{T, <:TriLower}, i::Integer, j::Integer, s::AbstractString) where T
	i >= j ? s : Base.replace_with_centered_mark(s)
end
function Base.replace_in_print_matrix(::TriMatrix{T, <:TriUpper}, i::Integer, j::Integer, s::AbstractString) where T
	i <= j ? s : Base.replace_with_centered_mark(s)
end


text:
using Base.Threads
@testset "rl_unplugged_dm" begin
    @testset "dm_control_suite" begin
        ds = rl_unplugged_dm_dataset(
            "fish_swim",
            [1, 2];
            type="dm_control_suite",
            is_shuffle = true,
            shuffle_buffer_size=10_000,
            tf_reader_bufsize=10_000,
            tf_reader_sz=10_000,
            batch_size=256,
            n_preallocations=nthreads()*12
        )

        @test typeof(ds)<:RingBuffer

        data = take!(ds)

        batch_size = 256
        feature_size = ReinforcementLearningDatasets.DM_CONTROL_SUITE_SIZE["fish_swim"]

        @test typeof(data.state) <: NamedTuple
        @test typeof(data.next_state) <: NamedTuple

        for feature in keys(feature_size)
            if split(feature, "/")[1] != "observation"
                if feature != "step_type"
                    ob_key = Symbol(feature)
                    @test size(getfield(data, ob_key)) == (feature_size[feature]..., batch_size,)
                end
            else
                state = data.state
                next_state = data.next_state
                ob_key = Symbol(chop(feature, head=length("observation")+1, tail=0))
                @test size(getfield(state, ob_key)) == (feature_size[feature]...,batch_size)
                @test size(getfield(next_state, ob_key)) == (feature_size[feature]..., batch_size,)
            end
        end
    end

    @testset "dm_locomotion_humanoid" begin
        ds = rl_unplugged_dm_dataset(
            "humanoid_corridor",
            [1, 2];
            type="dm_locomotion_humanoid",
            is_shuffle = true,
            shuffle_buffer_size=10_000,
            tf_reader_bufsize=10_000,
            tf_reader_sz=10_000,
            batch_size=256,
            n_preallocations=nthreads()*12
        )

        @test typeof(ds)<:RingBuffer

        data = take!(ds)

        batch_size = 256
        feature_size = ReinforcementLearningDatasets.DM_LOCOMOTION_HUMANOID_SIZE

        @test typeof(data.state) <: NamedTuple
        @test typeof(data.next_state) <: NamedTuple

        for feature in keys(feature_size)
            if split(feature, "/")[1] != "observation"
                if feature != "step_type"
                    ob_key = Symbol(feature)
                    @test size(getfield(data, ob_key)) == (feature_size[feature]..., batch_size,)
                end
            else
                state = data.state
                next_state = data.next_state
                ob_key = Symbol(chop(feature, head=length("observation")+1, tail=0))
                @test size(getfield(state, ob_key)) == (feature_size[feature]..., batch_size,)
                @test size(getfield(next_state, ob_key)) == (feature_size[feature]..., batch_size,)
            end
        end
    end

    @testset "dm_locomotion_rodent" begin
        ds = rl_unplugged_dm_dataset(
            "rodent_escape",
            [1, 2];
            type="dm_locomotion_rodent",
            is_shuffle = true,
            shuffle_buffer_size=10_000,
            tf_reader_bufsize=10_000,
            tf_reader_sz=10_000,
            batch_size=256,
            n_preallocations=nthreads()*12
        )

        @test typeof(ds)<:RingBuffer

        data = take!(ds)

        batch_size = 256
        feature_size = ReinforcementLearningDatasets.DM_LOCOMOTION_RODENT_SIZE

        @test typeof(data.state) <: NamedTuple
        @test typeof(data.next_state) <: NamedTuple

        for feature in keys(feature_size)
            if split(feature, "/")[1] != "observation"
                if feature != "step_type"
                    ob_key = Symbol(feature)
                    @test size(getfield(data, ob_key)) == (feature_size[feature]..., batch_size,)
                end
            else
                state = data.state
                next_state = data.next_state
                ob_key = Symbol(chop(feature, head=length("observation")+1, tail=0))
                @test size(getfield(state, ob_key)) == (feature_size[feature]..., batch_size,)
                @test size(getfield(next_state, ob_key)) == (feature_size[feature]..., batch_size,)
            end
        end
    end
end

ext: jl


text:
abstract type AbstractSelection end
(-)(s::Pair{<:AbstractSelection,S}) where S = -first(s) => last(s)
(!)(s::AbstractSelection) = -s
(~)(s::AbstractSelection) = -s
bool(s::AbstractSelection) = s.b
function Base.show(io::IO, s::AbstractSelection)
    !bool(s) && print(io, "-")
    print(io, typeof(s).name.name)
    :f in fieldnames(typeof(s)) ?
        print(io, "(", s.f, ")") :
        print(io, "(", s.s isa Symbol ? ":" : "", s.s, ")")
end

abstract type AbstractContextSelection <: AbstractSelection end
abstract type AbstractMultiSelection <: AbstractSelection end
Base.first(s::AbstractMultiSelection) = s.s1
Base.last(s::AbstractMultiSelection) = s.s2

abstract type AbstractMetaSelection end
abstract type AbstractRenaming end
# abstract type AbstractPermuteCols end
abstract type AbstractTransformation{Bool} end

function translate_types end
function selection end
function renaming end
function transformation end

function cols end
function not end


text:
abstract ClassificationModel

type MultiAdversarialModel <: ClassificationModel
  w::Vector{Float64}
  alpha::Vector{Float64}
  constraints::Vector{Tuple{Integer, Vector{Integer}}}
  n_class::Int
  game_value_01::Float64
  game_value_augmented::Float64
  train_adv_loss::Float64
  train_01_loss::Float64
end

type KernelMultiAdversarialModel <: ClassificationModel
  kernel::Symbol
  kernel_params::Vector{Float64}
  alpha::Vector{Float64}
  constraints::Vector{Tuple{Integer, Vector{Integer}}}
  n_class::Int
  game_value_01::Float64
  game_value_augmented::Float64
  train_adv_loss::Float64
  train_01_loss::Float64
end

max_stars_repo_name: Koprvhdix/adversarial-multiclass

max_issues_repo_name: Koprvhdix/adversarial-multiclass

max_forks_repo_name: Koprvhdix/adversarial-multiclass


text:
# MIT License

# Copyright (c) Microsoft Corporation.

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:

# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.

# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE

# OpticalInterface.jl is included in Geometry.jl so is omitted here

include("OpticalRay.jl")
include("Emitter.jl")
include("Emitter2.jl")
include("Fresnel.jl")
include("Paraxial.jl")
include("Grating.jl")
include("LensAssembly.jl")
include("HierarchicalImage.jl")
include("OpticalSystem.jl")
include("Lenses.jl")
include("Eye.jl")
include("Fields.jl")


text:
@testset "Normalized Cut" begin

    gx = SimpleGraph(6)
    add_edge!(gx, 1, 2)
    add_edge!(gx, 2, 3)
    add_edge!(gx, 1, 3)
    add_edge!(gx, 4, 5)
    add_edge!(gx, 4, 6)
    add_edge!(gx, 5, 6)
    add_edge!(gx, 1, 6)
    add_edge!(gx, 3, 4)

    w = zeros(6, 6)
    w[2, 1] = 1.0
    w[3, 1] = 1.0
    w[6, 1] = 0.1
    w[1, 2] = 1.0
    w[3, 2] = 1.0
    w[1, 3] = 1.0
    w[2, 3] = 1.0
    w[4, 3] = 0.2
    w[3, 4] = 0.2
    w[5, 4] = 1.0
    w[6, 4] = 1.0
    w[4, 5] = 1.0
    w[6, 5] = 1.0
    w[1, 6] = 0.1
    w[4, 6] = 1.0
    w[5, 6] = 1.0

    for g in testgraphs(gx)
      labels = @inferred(normalized_cut(g, 1, w))
      @test labels == [1, 1, 1, 2, 2, 2] || labels == [2, 2, 2, 1, 1, 1]
    end

    w = SparseMatrixCSC(w)
    for g in testgraphs(gx)
      labels = @inferred(normalized_cut(g, 1, w))
      @test labels == [1, 1, 1, 2, 2, 2] || labels == [2, 2, 2, 1, 1, 1]
    end

    gx = SimpleGraph(4)
    add_edge!(gx, 1, 2)
    add_edge!(gx, 3, 4)

    w = zeros(4, 4)
    w[2, 1] = 1.0
    w[1, 2] = 1.0
    w[4, 3] = 1.0
    w[3, 4] = 1.0
    for g in testgraphs(gx)
      labels = @inferred(normalized_cut(g, 0.1, w))
      @test labels == [1, 1, 2, 2] || labels == [2, 2, 1, 1]
    end

    w = SparseMatrixCSC(w)
    for g in testgraphs(gx)
      labels = @inferred(normalized_cut(g, 0.1, w))
      @test labels == [1, 1, 2, 2] || labels == [2, 2, 1, 1]
    end

    w = ones(12, 12)
    for g in testgraphs(gx)
      labels = @inferred(normalized_cut(g, 0.1, w))
      @test labels == [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    end

    w = ones(12, 12)
    for g in testgraphs(gx)
      labels = @inferred(normalized_cut(g, 0.1, w))
      @test labels == [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
    end

    g = PathGraph(50)

    function contiguous(labels::Vector{Int})::Bool
        changes = 0
        for i in 1:length(labels)-1
            if labels[i] != labels[i+1]
                changes += 1
            end
        end
        return changes == length(unique(labels)) - 1
    end

    num_subgraphs = Vector{Int}(undef, 9)

    for t in [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
        labels = @inferred(normalized_cut(g, t))
        @test contiguous(labels) == true
        num_subgraphs[convert(Int, 10*t)] = size(unique(labels), 1)
    end

    @test issorted(num_subgraphs) == true

    @test any(length(unique(@show normalized_cut(g,t))) == 4 for t in [0.125, 0.15, 0.16, 0.175, 0.20])

end


text:
module Models

using ..FastAI

using BSON
using Flux
using Zygote
using DataDeps
using Metalhead

include("./datadeps.jl")

function __init__()
    initdatadeps()
end

include("layers.jl")
include("blocks.jl")

include("./metalhead.jl")
include("./xresnet.jl")

export xresnet18, xresnet50

end


text:
"""
    AbstractInput

Abstract supertype for all input types.

### Notes

The input types defined here implement an iterator interface, such that other methods
can build upon the behavior of inputs which are either constant or varying.

Iteration is supported with an index number called *iterator state*.
The iteration function `Base.iterate` takes and returns a tuple (`input`, `state`),
where `input` represents the value of the input, and `state` is an index which
counts the number of times this iterator was called.

A convenience function `nextinput(input, n)` is also provided and it returns the
first `n` elements of `input`.
"""
abstract type AbstractInput end

"""
    ConstantInput{UT} <: AbstractInput

Type representing an input that remains constant in time.

### Fields

- `U` -- input set

### Examples

The constant input holds a single element and its length is infinite.
To access the field `U`, you can use Base's `iterate` given a state, or the method
 `nextinput` given the number of desired input elements:

```jldoctest constant_input
julia> c = ConstantInput(-1//2)
ConstantInput{Rational{Int64}}(-1//2)

julia> iterate(c, 1)
(-1//2, nothing)

julia> iterate(c, 2)
(-1//2, nothing)

julia> collect(nextinput(c, 4))
4-element Array{Rational{Int64},1}:
 -1//2
 -1//2
 -1//2
 -1//2
```

The elements of this input are rational numbers:

```jldoctest constant_input
julia> eltype(c)
Rational{Int64}
```

To transform a constant input, you can use `map` as in:

```jldoctest constant_input
julia> map(x->2*x, c)
ConstantInput{Rational{Int64}}(-1//1)
```
"""
struct ConstantInput{UT} <: AbstractInput
    U::UT
end

Base.eltype(::Type{ConstantInput{UT}}) where {UT} = UT
Base.iterate(input::ConstantInput, state::Union{Int, Nothing}=nothing) = (input.U, nothing)
Base.IteratorSize(::Type{<:ConstantInput}) = Base.IsInfinite()
Base.IteratorEltype(::Type{<:ConstantInput}) = Base.HasEltype()

Base.map(f::Function, c::ConstantInput) = ConstantInput(f(c.U))

"""
    nextinput(input::ConstantInput, n::Int=1)

Returns the first `n` elements of this input.

### Input

- `input` -- a constant input
- `n`     -- (optional, default: `1`) the number of desired elements

### Output

A repeated iterator that generates `n` equal samples of this input.
"""
nextinput(input::ConstantInput, n::Int=1) = Base.Iterators.repeated(input.U, n)

"""
    VaryingInput{UT, VUT<:AbstractVector{UT}} <: AbstractInput

Type representing an input that may vary with time.

### Fields

- `U` -- vector of input sets

### Examples

The varying input holds a vector and its length equals the number
of elements in the vector. Consider an input given by a vector of rational numbers:

```jldoctest varying_input
julia> v = VaryingInput([-1//2, 1//2])
VaryingInput{Rational{Int64},Array{Rational{Int64},1}}(Rational{Int64}[-1//2, 1//2])

julia> length(v)
2

julia> eltype(v)
Rational{Int64}
```

Base's `iterate` method receives the input and an integer state and returns the
input element and the next iteration state:

```jldoctest varying_input
julia> iterate(v, 1)
(-1//2, 2)

julia> iterate(v, 2)
(1//2, 3)
```

The method `nextinput` receives a varying input and an integer `n` and returns
an iterator over the first `n` elements of this input (where `n=1` by default):

```jldoctest varying_input
julia> typeof(nextinput(v))
Base.Iterators.Take{VaryingInput{Rational{Int64},Array{Rational{Int64},1}}}

julia> collect(nextinput(v, 1))
1-element Array{Rational{Int64},1}:
 -1//2

julia> collect(nextinput(v, 2))
2-element Array{Rational{Int64},1}:
 -1//2
  1//2
```

You can collect the inputs in an array, or equivalently use list comprehension,
(or use a `for` loop):

```jldoctest varying_input
julia> collect(v)
2-element Array{Rational{Int64},1}:
 -1//2
  1//2

julia> [2*vi for vi in v]
2-element Array{Rational{Int64},1}:
 -1//1
  1//1
```

Since this input type is finite, querying more elements than its length returns
the full vector:

```jldoctest varying_input
julia> collect(nextinput(v, 4))
2-element Array{Rational{Int64},1}:
 -1//2
  1//2
```

To transform a varying input, you can use `map` as in:

```jldoctest varying_input
julia> map(x->2*x, v)
VaryingInput{Rational{Int64},Array{Rational{Int64},1}}(Rational{Int64}[-1//1, 1//1])
```
"""
struct VaryingInput{UT, VUT<:AbstractVector{UT}} <: AbstractInput
    U::VUT  # input sequence
end

Base.eltype(::Type{VaryingInput{UT, VUT}}) where {UT, VUT} = UT

function Base.iterate(input::VaryingInput, state::Int=1)
    if state > length(input.U)
        return nothing
    else
        return (input.U[state], state + 1)
    end
end

Base.length(input::VaryingInput) = length(input.U)

Base.IteratorSize(::Type{<:VaryingInput}) = Base.HasLength()
Base.IteratorEltype(::Type{<:VaryingInput}) = Base.HasEltype()

Base.map(f::Function, v::VaryingInput) = VaryingInput(f.(v.U))

"""
    nextinput(input::VaryingInput, n::Int=1)

Returns the first `n` elements of this input.

### Input

- `input` -- varying input
- `n`     -- (optional, default: `1`) number of desired elements

### Output

An iterator of type `Base.Iterators.Take` that represents at most the first
`n` elements of this input.
"""
nextinput(input::VaryingInput, n::Int=1) = Base.Iterators.take(input, n)


text:
@testset "Parallel.Betweenness" begin

    s2 = SimpleDiGraph(3)
    add_edge!(s2, 1, 2); add_edge!(s2, 2, 3); add_edge!(s2, 3, 3)
    s1 = SimpleGraph(s2)
    g3 = PathGraph(5)
    gint = loadgraph(joinpath(testdir, "testdata", "graph-50-500.jgz"), "graph-50-500")
    for g in testdigraphs(gint)
        z  = @inferred(LightGraphs.betweenness_centrality(g))
        zt = @inferred(Parallel.betweenness_centrality(g; parallel=:threads))
        @test all(isapprox.(z, zt))
        zd = @inferred(Parallel.betweenness_centrality(g; parallel=:distributed))
        @test all(isapprox.(z, zd))

        y  = LightGraphs.betweenness_centrality(g, endpoints=true, normalize=false)
        yt = Parallel.betweenness_centrality(g, endpoints=true, normalize=false, parallel=:threads)
        @test all(isapprox(y, yt))
        yd = Parallel.betweenness_centrality(g, endpoints=true, normalize=false, parallel=:distributed)
        @test all(isapprox(y, yd))

        xt = @inferred(Parallel.betweenness_centrality(g, 3; parallel=:threads))
        @test length(xt) == 50
        xd = @inferred(Parallel.betweenness_centrality(g, 3; parallel=:distributed))
        @test length(xd) == 50

        xt2 = @inferred(Parallel.betweenness_centrality(g, collect(1:20); parallel=:threads))
        @test length(xt2) == 50
        xd2 = @inferred(Parallel.betweenness_centrality(g, collect(1:20); parallel=:distributed))
        @test length(xd2) == 50
    end
    @test @inferred(Parallel.betweenness_centrality(s1; parallel=:threads)) == [0, 1, 0]
    @test @inferred(Parallel.betweenness_centrality(s1; parallel=:distributed)) == [0, 1, 0]

    @test @inferred(Parallel.betweenness_centrality(s2; parallel=:threads)) == [0, 0.5, 0]
    @test @inferred(Parallel.betweenness_centrality(s2; parallel=:distributed)) == [0, 0.5, 0]

    g = SimpleGraph(2)
    add_edge!(g, 1, 2)
    z  = LightGraphs.betweenness_centrality(g; normalize=true)
    @test z[1] == z[2] == 0.0
    zt = Parallel.betweenness_centrality(g; normalize=true, parallel=:threads)
    @test all(isapprox(z, zt))
    zd = Parallel.betweenness_centrality(g; normalize=true, parallel=:distributed)
    @test all(isapprox(z, zd))

    z2  = LightGraphs.betweenness_centrality(g, vertices(g))
    zt2 = Parallel.betweenness_centrality(g, vertices(g); parallel=:threads)
    @test all(isapprox(z2, zt2))
    zd2 = Parallel.betweenness_centrality(g, vertices(g); parallel=:distributed)
    @test all(isapprox(z2, zd2))

    z3  = LightGraphs.betweenness_centrality(g, [vertices(g);])
    zt3 = Parallel.betweenness_centrality(g, [vertices(g);]; parallel=:threads)
    @test all(isapprox(z3, zt3))
    zd3 = Parallel.betweenness_centrality(g, [vertices(g);]; parallel=:distributed)
    @test all(isapprox(z3, zd3))

    @test z == z2 == z3

    z  = LightGraphs.betweenness_centrality(g3; normalize=false)
    zt = Parallel.betweenness_centrality(g3; normalize=false, parallel=:threads)
    @test all(isapprox(z, zt))
    zd = Parallel.betweenness_centrality(g3; normalize=false, parallel=:distributed)
    @test all(isapprox(z, zd))
    @test z[1] == z[5] == 0.0
end


text:
##############
# Reflection #
##############

mutable struct Reflection
    signature::DataType
    method::Method
    static_params::Vector{Any}
    code_info::CodeInfo
end

# Return `Reflection` for signature `sigtypes` and `world`, if possible. Otherwise, return `nothing`.
function reflect(@nospecialize(sigtypes::Tuple), world::UInt = typemax(UInt))
    if length(sigtypes) > 2 && sigtypes[1] === typeof(invoke)
        @assert sigtypes[3] <: Type{<:Tuple}
        sigtypes = (sigtypes[2], sigtypes[3].parameters[1].parameters...)
    end
    # This works around a subtyping bug. Basically, callers can deconstruct upstream
    # `UnionAll` types in such a way that results in a type with free type variables, in
    # which case subtyping can just break.
    #
    # God help you if you try to use a type parameter here (e.g. `::Type{S} where S<:Tuple`)
    # instead of this nutty workaround, because the compiler can just rewrite `S` into
    # whatever it thinks is "type equal" to the actual provided value. In other words, if
    # `S` is defined as e.g. `f(::Type{S}) where S`, and you call `f(T)`, you should NOT
    # assume that `S === T`. If you did, SHAME ON YOU. It doesn't matter that such an
    # assumption holds true for essentially all other kinds of values. I haven't counted in
    # a while, but I'm pretty sure I have ~40+ hellish years of Julia experience, and this
    # still catches me every time. Who even uses this crazy language?
    S = Tuple{map(s -> Core.Compiler.has_free_typevars(s) ? typeof(s.parameters[1]) : s, sigtypes)...}
    (S.parameters[1]::DataType).name.module === Core.Compiler && return nothing
    _methods = Base._methods_by_ftype(S, -1, world)
    method_index = 0
    for i in 1:length(_methods)
        if _methods[i][1] === S
            method_index = i
            break
        end
    end
    method_index === 0 && return nothing
    type_signature, raw_static_params, method = _methods[method_index]
    method_instance = Core.Compiler.code_for_method(method, type_signature, raw_static_params, world, false)
    method_instance === nothing && return nothing
    method_signature = method.sig
    static_params = Any[raw_static_params...]
    code_info = Core.Compiler.retrieve_code_info(method_instance)
    isa(code_info, CodeInfo) || return nothing
    code_info = Core.Compiler.copy_code_info(code_info)
    return Reflection(S, method, static_params, code_info)
end

###########
# overdub #
###########

const OVERDUB_CTX_SYMBOL = gensym("overdub_context")
const OVERDUB_ARGS_SYMBOL = gensym("overdub_arguments")
const OVERDUB_TMP_SYMBOL = gensym("overdub_tmp")

# The `overdub` pass has four intertwined tasks:
#   1. Apply the user-provided pass, if one is given
#   2. Munge the reflection-generated IR into a valid form for returning from
#      `recurse_generator` (i.e. add new argument slots, substitute static
#      parameters, destructure overdub arguments into underlying method slots, etc.)
#   3. Replace all calls of the form `output = f(args...)` with:
#      ```
#      prehook(ctx, f, args...)
#      tmp = execute(ctx, f, args...)
#      isa(tmp, OverdubInstead) ? overdub(ctx, f, args...) : tmp
#      posthook(ctx, f, args...)
#      output = tmp
#      ```
#   4. If tagging is enabled, do the necessary IR transforms for the metadata tagging system
function overdub_pass!(reflection::Reflection,
                       context_type::DataType,
                       pass_type::DataType = NoPass,
                       is_invoke::Bool = false)
    signature = reflection.signature
    method = reflection.method
    static_params = reflection.static_params
    code_info = reflection.code_info

    # TODO: This `iskwfunc` is part of a hack that `overdub_pass!` implements in order to fix
    # jrevels/Cassette.jl#48. The assumptions made by this hack are quite fragile, so we
    # should eventually get Base to expose a standard/documented API for this. Here, we see
    # this hack's first assumption: that `Core.kwfunc(f)` is going to return a function whose
    # type name is prefixed by `#kw##`. More assumptions for this hack will be commented on
    # as we go.
    iskwfunc = startswith(String(signature.parameters[1].name.name), "#kw##")
    istaggingenabled = hastagging(context_type)

    #=== execute user-provided pass (is a no-op by default) ===#

    if !iskwfunc
        code_info = pass_type(context_type, signature, code_info)
    end

    #=== munge the code into a valid form for `overdub_generator` ===#

    # NOTE: The slotflags set by this pass are set according to what makes sense based on the
    # compiler's actual `@code_lowered` output in practice, since this real-world output does
    # not seem to match Julia's developer documentation.

    # construct new slotnames/slotflags for added slots
    code_info.slotnames = Any[:overdub, OVERDUB_CTX_SYMBOL, OVERDUB_ARGS_SYMBOL, code_info.slotnames..., OVERDUB_TMP_SYMBOL]
    code_info.slotflags = UInt8[0x00, 0x00, 0x00, code_info.slotflags..., 0x00]
    n_prepended_slots = 3
    overdub_ctx_slot = SlotNumber(2)
    overdub_args_slot = SlotNumber(3)
    overdub_tmp_slot = SlotNumber(length(code_info.slotnames))
    invoke_offset = is_invoke ? 2 : 0

    # For the sake of convenience, the rest of this pass will translate `code_info`'s fields
    # into these overdubbed equivalents instead of updating `code_info` in-place. Then, at
    # the end of the pass, we'll reset `code_info` fields accordingly.
    overdubbed_code = Any[]
    overdubbed_codelocs = Int32[]

    # destructure the generated argument slots into the overdubbed method's argument slots.
    n_actual_args = fieldcount(signature)
    n_method_args = Int(method.nargs)
    for i in 1:n_method_args
        slot = i + n_prepended_slots
        actual_argument = Expr(:call, GlobalRef(Core, :getfield), overdub_args_slot, i + invoke_offset)
        push!(overdubbed_code, :($(SlotNumber(slot)) = $actual_argument))
        push!(overdubbed_codelocs, code_info.codelocs[1])
        code_info.slotflags[slot] |= 0x02 # ensure this slotflag has the "assigned" bit set
    end

    # If `method` is a varargs method, we have to restructure the original method call's
    # trailing arguments into a tuple and assign that tuple to the expected argument slot.
    if method.isva
        if !isempty(overdubbed_code)
            # remove the final slot reassignment leftover from the previous destructuring
            pop!(overdubbed_code)
            pop!(overdubbed_codelocs)
        end
        if hastagging(context_type)
            trailing_arguments = Expr(:call, GlobalRef(Cassette, :_tagged_new_tuple_unsafe), overdub_ctx_slot)
        else
            trailing_arguments = Expr(:call, GlobalRef(Core, :tuple))
        end
        for i in n_method_args:n_actual_args
            push!(overdubbed_code, Expr(:call, GlobalRef(Core, :getfield), overdub_args_slot, i + invoke_offset))
            push!(overdubbed_codelocs, code_info.codelocs[1])
            push!(trailing_arguments.args, SSAValue(length(overdubbed_code)))
        end
        push!(overdubbed_code, Expr(:(=), SlotNumber(n_method_args + n_prepended_slots), trailing_arguments))
        push!(overdubbed_codelocs, code_info.codelocs[1])
    end

    #=== finish initialization of `overdubbed_code`/`overdubbed_codelocs` ===#

    # substitute static parameters, offset slot numbers by number of added slots, and
    # offset statement indices by the number of additional statements
    Base.Meta.partially_inline!(code_info.code, Any[], method.sig, static_params,
                                n_prepended_slots, length(overdubbed_code), :propagate)

    original_code_start_index = length(overdubbed_code) + 1

    append!(overdubbed_code, code_info.code)
    append!(overdubbed_codelocs, code_info.codelocs)

    #=== perform tagged module transformation if tagging is enabled ===#

    if istaggingenabled && !iskwfunc
        # find `GlobalRef`s in IR and set up tagged module replacements
        modules = Any[]
        original_code_region = view(overdubbed_code, original_code_start_index:length(overdubbed_code))
        replace_match!(x -> isa(x, GlobalRef), original_code_region) do x
            m = GlobalRef(parentmodule(x.mod), nameof(x.mod))
            i = findfirst(isequal(m), modules)
            if isa(i, Nothing)
                push!(modules, m)
                i = length(modules)
            end
            return Expr(:replaceglobalref, original_code_start_index + i - 1, x)
        end
        for i in 1:length(modules)
            modules[i] = Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :fetch_tagged_module)), overdub_ctx_slot, modules[i])
        end

        # insert `fetch_tagged_module`s at the `original_code_start_index`
        insert_statements!(overdubbed_code, overdubbed_codelocs,
                            (x, i) -> i == original_code_start_index ? length(modules) + 1 : nothing,
                            (x, i) -> [modules..., x])

        # append `tagged_globalref_set_meta!` to `GlobalRef() = ...` statements
        insert_statements!(overdubbed_code, overdubbed_codelocs,
                            (x, i) -> begin
                                if (i > original_code_start_index &&
                                    Base.Meta.isexpr(x, :(=)) &&
                                    Base.Meta.isexpr(x.args[1], :replaceglobalref))
                                    return 3
                                end
                                return nothing
                            end,
                            (x, i) -> begin
                                lhs, rhs = x.args
                                tagmodssa = SSAValue(lhs.args[1])
                                globalref = lhs.args[2]
                                name = QuoteNode(globalref.name)
                                return [
                                    rhs,
                                    Expr(:(=), globalref, Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :untag)), SSAValue(i), overdub_ctx_slot)),
                                    Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :tagged_globalref_set_meta!)), overdub_ctx_slot, tagmodssa, name, SSAValue(i))
                                ]
                            end)

        # replace `GlobalRef`-loads with `Cassette.tagged_globalref`
        original_code_start_index += length(modules)
        insert_statements!(overdubbed_code, overdubbed_codelocs,
                            (x, i) -> begin
                                i >= original_code_start_index || return nothing
                                stmt = Base.Meta.isexpr(x, :(=)) ? x.args[2] : x
                                Base.Meta.isexpr(stmt, :replaceglobalref) && return 1
                                if isa(stmt, Expr) # Base.Meta.isexpr(stmt, :call) || Base.Meta.isexpr(stmt, :new) || Base.Meta.isexpr(stmt, :return)
                                    count = 0
                                    for arg in stmt.args
                                        if Base.Meta.isexpr(arg, :replaceglobalref)
                                            count += 1
                                        end
                                    end
                                    count > 0 && return count + 1
                                end
                                return nothing
                            end,
                            (x, i) -> begin
                                items = Any[]
                                stmt = Base.Meta.isexpr(x, :(=)) ? x.args[2] : x
                                if Base.Meta.isexpr(stmt, :replaceglobalref)
                                    tagmodssa = SSAValue(stmt.args[1])
                                    globalref = stmt.args[2]
                                    name = QuoteNode(globalref.name)
                                    result = Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :tagged_globalref)), overdub_ctx_slot, tagmodssa, name, globalref)
                                elseif isa(stmt, Expr) # Base.Meta.isexpr(stmt, :call) || Base.Meta.isexpr(stmt, :new) || Base.Meta.isexpr(stmt, :return)
                                    result = Expr(stmt.head)
                                    for arg in stmt.args
                                        if Base.Meta.isexpr(arg, :replaceglobalref)
                                            tagmodssa = SSAValue(arg.args[1])
                                            globalref = arg.args[2]
                                            name = QuoteNode(globalref.name)
                                            push!(result.args, SSAValue(i + length(items)))
                                            push!(items, Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :tagged_globalref)), overdub_ctx_slot, tagmodssa, name, globalref))
                                        else
                                            push!(result.args, arg)
                                        end
                                    end
                                end
                                if Base.Meta.isexpr(x, :(=))
                                    result = Expr(:(=), x.args[1], result)
                                end
                                push!(items, result)
                                return items
                            end)
    end

    #=== untag all `foreigncall` SSAValue/SlotNumber arguments if tagging is enabled ===#

    if istaggingenabled && !iskwfunc
        insert_statements!(overdubbed_code, overdubbed_codelocs,
                            (x, i) -> begin
                                stmt = Base.Meta.isexpr(x, :(=)) ? x.args[2] : x
                                if Base.Meta.isexpr(stmt, :foreigncall)
                                    count = 0
                                    for arg in stmt.args
                                        if isa(arg, SSAValue) || isa(arg, SlotNumber)
                                            count += 1
                                        end
                                    end
                                    count > 0 && return count + 1
                                end
                                return nothing
                            end,
                            (x, i) -> begin
                                items = Any[]
                                stmt = Base.Meta.isexpr(x, :(=)) ? x.args[2] : x
                                result = Expr(:foreigncall)
                                for arg in stmt.args
                                    if isa(arg, SSAValue) || isa(arg, SlotNumber)
                                        push!(result.args, SSAValue(i + length(items)))
                                        push!(items, Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :untag)), arg, overdub_ctx_slot))
                                    else
                                        push!(result.args, arg)
                                    end
                                end
                                if Base.Meta.isexpr(x, :(=))
                                    result = Expr(:(=), x.args[1], result)
                                end
                                push!(items, result)
                                return items
                            end)
    end

    #=== untag `gotoifnot` conditionals if tagging is enabled ===#

    if istaggingenabled && !iskwfunc
        insert_statements!(overdubbed_code, overdubbed_codelocs,
                            (x, i) -> Base.Meta.isexpr(x, :gotoifnot) ? 2 : nothing,
                            (x, i) -> [
                                Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :untag)), x.args[1], overdub_ctx_slot),
                                Expr(:gotoifnot, SSAValue(i), x.args[2])
                            ])
    end

    #=== replace `Expr(:new, ...)` with `Expr(:call, :tagged_new)` if tagging is enabled ===#

    if istaggingenabled && !iskwfunc
        replace_match!(x -> Base.Meta.isexpr(x, :new), overdubbed_code) do x
            return Expr(:call, Expr(:nooverdub, GlobalRef(Cassette, :tagged_new)), overdub_ctx_slot, x.args...)
        end
    end

    #=== replace `Expr(:call, ...)` with `Expr(:call, :overdub, ...)` calls ===#

    if iskwfunc
        # Another assumption of this `iskwfunc` hack is that the second to last statement in
        # the lowered IR for `Core.kwfunc(f)` is the call to the "underlying" non-kwargs form
        # of `f`. Thus, we `overdub` that call instead of replacing it with `call`.
        for i in 1:length(overdubbed_code)
            stmt = overdubbed_code[i]
            replacewith = i === (length(overdubbed_code) - 1) ? :overdub : :call
            if Base.Meta.isexpr(stmt, :(=))
                replacein = stmt.args
                replaceat = 2
            else
                replacein = overdubbed_code
                replaceat = i
            end
            stmt = replacein[replaceat]
            if Base.Meta.isexpr(stmt, :call)
                replacein[replaceat] = Expr(:call, GlobalRef(Cassette, replacewith), overdub_ctx_slot, stmt.args...)
            end
        end
    else
        stmtcount = (x, i) -> begin
            i >= original_code_start_index || return nothing
            stmt = Base.Meta.isexpr(x, :(=)) ? x.args[2] : x
            if Base.Meta.isexpr(stmt, :call) && !(Base.Meta.isexpr(stmt.args[1], :nooverdub))
                return 7
            end
            return nothing
        end
        newstmts = (x, i) -> begin
            callstmt = Base.Meta.isexpr(x, :(=)) ? x.args[2] : x
            execstmt = Expr(:call, GlobalRef(Cassette, :execute), overdub_ctx_slot, callstmt.args...)
            overdubstmt = Expr(:call, GlobalRef(Cassette, :overdub), overdub_ctx_slot, callstmt.args...)
            return [
                Expr(:call, GlobalRef(Cassette, :prehook), overdub_ctx_slot, callstmt.args...),
                Expr(:(=), overdub_tmp_slot, execstmt),
                Expr(:call, GlobalRef(Core, :isa), overdub_tmp_slot, GlobalRef(Cassette, :OverdubInstead)),
                Expr(:gotoifnot, SSAValue(i + 2), i + 5),
                Expr(:(=), overdub_tmp_slot, overdubstmt),
                Expr(:call, GlobalRef(Cassette, :posthook), overdub_ctx_slot, overdub_tmp_slot, callstmt.args...),
                Base.Meta.isexpr(x, :(=)) ? Expr(:(=), x.args[1], overdub_tmp_slot) : overdub_tmp_slot
            ]
        end
        insert_statements!(overdubbed_code, overdubbed_codelocs, stmtcount, newstmts)
    end

    #=== unwrap all `Expr(:nooverdub)`s ===#

    replace_match!(x -> x.args[1], x -> Base.Meta.isexpr(x, :nooverdub), overdubbed_code)

    #=== replace all `Expr(:contextslot)`s ===#

    replace_match!(x -> overdub_ctx_slot, x -> Base.Meta.isexpr(x, :contextslot), overdubbed_code)

    #=== set `code_info`/`reflection` fields accordingly ===#

    code_info.code = overdubbed_code
    code_info.codelocs = overdubbed_codelocs
    code_info.ssavaluetypes = length(overdubbed_code)
    code_info.method_for_inference_limit_heuristics = method
    reflection.code_info = code_info

    return reflection
end

# `args` is `(typeof(original_function), map(typeof, original_args_tuple)...)`
function __overdub_generator__(pass_type, self, context_type, args::Tuple)
    if nfields(args) > 0
        is_builtin = args[1] <: Core.Builtin
        is_invoke = args[1] === typeof(Core.invoke)
        if !is_builtin || is_invoke
            try
                untagged_args = ((untagtype(args[i], context_type) for i in 1:nfields(args))...,)
                reflection = reflect(untagged_args)
                if isa(reflection, Reflection)
                    overdub_pass!(reflection, context_type, pass_type, is_invoke)
                    return reflection.code_info
                end
            catch err
                errmsg = "ERROR COMPILING $args IN CONTEXT $(context_type): \n" * sprint(showerror, err)
                return quote
                    error($errmsg)
                end
            end
        end
    end
    return quote
        $(Expr(:meta, :inline))
        $Cassette.fallback($OVERDUB_CTX_SYMBOL, $OVERDUB_ARGS_SYMBOL...)
    end
end

@inline apply_args(::ContextWithTag{Nothing}, args...) = Core._apply(Core.tuple, args...)
@inline apply_args(ctx::Context, args...) = tagged_apply_args(ctx, args...)

function overdub end

function overdub_definition(pass, line, file)
    return quote
        function $Cassette.overdub($OVERDUB_CTX_SYMBOL::$Cassette.ContextWithPass{pass}, $OVERDUB_ARGS_SYMBOL...) where {pass<:$pass}
            $(Expr(:meta,
                   :generated,
                   Expr(:new,
                        Core.GeneratedFunctionStub,
                        :__overdub_generator__,
                        Any[:overdub, OVERDUB_CTX_SYMBOL, OVERDUB_ARGS_SYMBOL],
                        Any[:pass],
                        line,
                        QuoteNode(Symbol(file)),
                        true)))
        end
        @inline function $Cassette.overdub(ctx::$Cassette.ContextWithPass{pass}, ::typeof(Core._apply), f, _args...) where {pass<:$pass}
            args = $Cassette.apply_args(ctx, _args...)
            $Cassette.prehook(ctx, f, args...)
            output = $Cassette.execute(ctx, f, args...)
            output = isa(output, $Cassette.OverdubInstead) ? $Cassette.overdub(ctx, f, args...) : output
            $Cassette.posthook(ctx, output, f, args...)
            return output
        end
    end
end

@eval $(overdub_definition(:NoPass, @__LINE__, @__FILE__))

@doc(
"""
```
overdub(context::Context, f, args...)
```

Execute `f(args...)` overdubbed with respect to `context`.

More specifically, execute `f(args...)`, but with every internal method invocation `g(x...)`
replaced by statements similar to the following:

```
begin
    prehook(context, g, x...)
    tmp = execute(context, g, x...)
    tmp = isa(tmp, Cassette.OverdubInstead) ? overdub(context, g, x...) : tmp
    posthook(context, tmp, g, x...)
    tmp
end
```

If Cassette cannot retrieve lowered IR for the method body of `f(args...)` (as determined by
`canoverdub(context, f, args...)`), then `overdub(context, f, args...)` will directly
translate to a call to `fallback(context, f, args...)`.

Additionally, for every method body encountered in execute trace, apply the compiler pass
associated with `context` if one exists. Note that this user-provided pass is performed on
the method IR before method invocations are transformed into the form specified above. See
the [`@pass`](@ref) macro for further details.

If `Cassette.hastagging(typeof(context))`, then a number of additional passes are run in
order to accomodate tagged value propagation:

- `Expr(:new)` is replaced with a call to `Cassette.tagged_new`
- conditional values passed to `Expr(:gotoifnot)` are untagged
- arguments to `Expr(:foreigncall)` are untagged
- load/stores to external module bindings are intercepted by the tagging system
""",
overdub)

"""
```
Cassette.@overdub(ctx, expression)
```

A convenience macro for executing `expression` within the context `ctx`. This macro roughly
expands to `Cassette.overdub(ctx, () -> expression)`.

See also: [`overdub`](@ref)
"""
macro overdub(ctx, expr)
    return :($Cassette.overdub($(esc(ctx)), () -> $(esc(expr))))
end


text:
const datafile = joinpath(@__DIR__, "..", "data", "DICE2010_082710d.xlsx")

"""
    read_params(f, range::String, count::Int, sheet::String="Base")

Get parameters from DICE2010 excel sheet.

`range` is a single cell or a range of cells in the excel sheet.
  Must be a cell reference of the form "A27" or a range "B56:B77".

`count` is the length of the time dimension; ignored if range
   refers to a single cell.

`sheet` is the name of the worksheet in the Excel file to read from.
  Defaults to "Base".

Examples:
    values = read_params(f, "B15:BI15", 40)   # read only the first 40 values

    value = read_params(f, "A27", sheet="Parameters")
    value = read_params(f, "A27:A27", sheet="Parameters") # same as above
"""
function read_params(f, range::String, T::Int=60; sheet::String="Base")
    data = f[sheet][range]
    parts = split(range, ":")
    return (length(parts) == 1 || parts[1] == parts[2]) ? data : Vector{Float64}(data[1:T])
end

function dice2010_excel_parameters(filename=datafile; nsteps=nothing)
    p = Dict{Symbol,Any}()

    # the number of time-indexed values to read
    nsteps = nsteps == nothing ? length(2005:10:2595) : nsteps

    # Open DICE_2010 Excel file to read parameters
    f = readxlsx(filename)

    # Preferences
    p[:elasmu]  = read_params(f, "B19")               # Elasticity of marginal utility of consumption
    p[:rr]      = read_params(f, "B18:BI18", nsteps)       # Social time preference factor

    # Population and technology
    p[:gama]    = read_params(f, "B9")                # Capital elasticity in production function
    p[:l]       = read_params(f, "B27:BI27", nsteps)       # Population (millions)
    p[:dk]      = read_params(f, "B10")               # Depreciation rate on capital (per year)
    p[:k0]      = read_params(f, "B13")               # Initial capital value (trill 2005 USD)
    p[:al]      = read_params(f, "B21:BI21", nsteps)       # Level of total factor productivity

    # Emissions parameters
    p[:sigma]   = read_params(f, "B46:BI46", nsteps)
    p[:etree]   = read_params(f, "B52:BI52", nsteps)
    p[:MIU] = read_params(f, "B133:BI133", nsteps)         # emissions control rate
    p[:S] = read_params(f, "B132:BI132", nsteps) / 100   # savings rate

    # Carbon cycle
    p[:mat0]    = read_params(f, "B57")   # Initial Concentration in atmosphere 2000 (GtC)
    p[:mat1]    = read_params(f, "B58")   # Initial Concentration in atmosphere 2010 (GtC)
    p[:mu0]     = read_params(f, "B59")   # Initial Concentration in biosphere/shallow oceans (GtC)
    p[:ml0]     = read_params(f, "B60")   # Initial Concentration in deep oceans (GtC)

    # Flow parameters
    p[:b12] = read_params(f, "B64") / 100 # Carbon cycle transition matrix atmosphere to shallow ocean
    p[:b23] = read_params(f, "B67") / 100 # Carbon cycle transition matrix shallow to deep ocean

    # Parameters for long-run consistency of carbon cycle
    p[:b11] = read_params(f, "B62") / 100 # Carbon cycle transition matrix atmosphere to atmosphere
    p[:b21] = read_params(f, "B63") / 100 # Carbon cycle transition matrix biosphere/shallow oceans to atmosphere
    p[:b22] = read_params(f, "B65") / 100 # Carbon cycle transition matrix shallow ocean to shallow oceans
    p[:b32] = read_params(f, "B66") / 100 # Carbon cycle transition matrix deep ocean to shallow ocean
    p[:b33] = read_params(f, "B68") / 100 # Carbon cycle transition matrix deep ocean to deep oceans

    # Climate model parameters
    p[:t2xco2]  = read_params(f, "B76")   # Equilibrium temp impact (oC per doubling CO2)
    p[:tatm0]   = read_params(f, "B73")   # Initial lower stratum temp change (C from 1900)
    p[:tatm1]   = read_params(f, "C73")   # Initial atmospheric temp change 2015 (C from 1900)
    p[:tocean0] = read_params(f, "B74")   # Initial lower stratum temp change (C from 1900)

    # Transient TSC Correction ("Speed of Adjustment Parameter")
    p[:c1]      = read_params(f, "B75")   # Climate equation coefficient for upper level
    p[:c3]      = read_params(f, "B78")   # Transfer coefficient upper to lower stratum
    p[:c4]      = read_params(f, "B79")   # Transfer coefficient for lower level
    p[:fco22x]  = read_params(f, "B77")   # Forcings of equilibrium CO2 doubling (Wm-2)

    # Climate damage parameters
    p[:a1] = read_params(f, "B33")   # Damage coefficient
    p[:a2] = read_params(f, "B34")   # Damage quadratic term
    p[:a3] = read_params(f, "B35")   # Damage exponent

    # Abatement cost
    p[:expcost2]    = read_params(f, "B44")               # Exponent of control cost function
    p[:pbacktime]   = read_params(f, "B42:BI42", nsteps)       # backstop price
    # Adjusted cost for backstop (or: "Abatement cost function coefficient")
    p[:cost1]       = read_params(f, "B37:BI37", nsteps)

    # Scaling and inessential parameters
    p[:scale1] = read_params(f, "B88")    # Multiplicative scaling coefficient
    p[:scale2] = read_params(f, "B89")    # Additive scaling coefficient

    # Exogenous forcing for other greenhouse gases
    p[:forcoth] = read_params(f, "B70:BI70", nsteps)

    # Fraction of emissions in control regime
    p[:partfract] = read_params(f, "B82:BI82", nsteps)

    # SLR Parameters

    p[:b1]    = read_params(f, "B51", sheet="Parameters")
    p[:b2]  = read_params(f, "B52", sheet="Parameters")
    p[:b3]      = read_params(f, "B53", sheet="Parameters")

    p[:therm0]      = read_params(f, "B178") # meters above 2000
    p[:gsic0]       = read_params(f, "B179")
    p[:gis0]        = read_params(f, "B180")
    p[:ais0]        = read_params(f, "B181")

    p[:therm_asym]  = read_params(f, "B173")
    p[:gsic_asym]   = read_params(f, "B174")
    p[:gis_asym]    = read_params(f, "B175")
    p[:ais_asym]    = read_params(f, "B176")

    p[:thermrate]   = read_params(f, "C173")
    p[:gsicrate]    = read_params(f, "C174")
    p[:gisrate]     = read_params(f, "C175")
    p[:aisrate]     = read_params(f, "C176")

    p[:slrthreshold] = read_params(f, "D176")

    return p
end


text:
using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using Test
using LinearAlgebra

@testset "GUN (native)" begin
    nep = nep_gallery("nlevp_native_gun")
    n = size(nep, 1)
    tol = 1e-11
    ref_eigenvalue = 22345.116783765 + 0.644998598im # from NLEIGS

    @bench @testset "Running algorithm" begin
        λ1,v1 = quasinewton(nep, λ = 150^2+1im, v = ones(n), logger = displaylevel, tol = tol, maxit = 500)

        v1 = v1 / norm(v1)

        @test norm(compute_Mlincomb(nep, λ1, v1)) < sqrt(tol)
        @test norm(compute_Mder(nep, λ1) * v1) < sqrt(tol)

        @test norm(λ1 - ref_eigenvalue) < sqrt(tol)*100
    end

    @bench @testset "Compute derivatives" begin
        λ = 150^2+2im
        v = randn(n)
        z1 = compute_Mlincomb(nep, λ, v, [1.0], 1)

        # Compare with divided difference
        ee = 1e-4
        z2 = (compute_Mlincomb(nep, λ+ee, v) - compute_Mlincomb(nep, λ-ee, v)) / (2*ee)

        @test norm(z2-z1) < (ee^2)*1000
    end

    @testset "Errors thrown" begin
        @test_throws MethodError nep_gallery("nlevp_native_gun", 15)
        @test_throws MethodError nep_gallery("nlevp_native_gun", t=15)
    end
end


text:
@block PietroVertechi [statsmakie] begin
    @cell "StatsMakie" [Group, Style, scatter, plot, wireframe, statistics] begin
        using StatsMakie
        using DataFrames, RDatasets
        using StatsMakie: linear, smooth
        N = 1000
        a = rand(1:2, N) # a discrete variable
        b = rand(1:2, N) # a discrete variable
        x = randn(N) # a continuous variable
        y = @. x * a + 0.8*randn() # a continuous variable
        z = x .+ y # a continuous variable
        @substep

        scatter(x, y, markersize = 0.2)
        @substep

        scatter(Group(a), x, y, markersize = 0.2)
        @substep

        scatter(Group(a), x, y, color = [:black, :red], markersize = 0.2)
        @substep

        scatter(Group(marker = a), x, y, markersize = 0.2)
        @substep

        scatter(Group(marker = a, color = b), x, y, markersize = 0.2)
        @substep

        scatter(Group(marker = a), Style(color = z), x, y)
        @substep

        scatter(Group(color = a), x, y, Style(markersize = z ./ 10))
        @substep

        plot(linear, x, y)
        @substep

        plot(linear, Group(a), x, y)
        @substep

        scatter(Group(a), x, y, markersize = 0.2)
        plot!(linear, Group(a), x, y)
        @substep

        plot(linear, Group(linestyle = a), x, y)
        @substep

        N = 200
        x = 10 .* rand(N)
        a = rand(1:2, N)
        y = sin.(x) .+ 0.5 .* rand(N) .+ cos.(x) .* a
        @substep

        scatter(Group(a), x, y)
        plot!(smooth, Group(a), x, y)
        @substep

        plot(histogram, y)
        @substep

        plot(histogram, x, y)

        @substep

        plot(histogram(nbins = 30), x, y)
        @substep

        wireframe(histogram(nbins = 30), x, y)
        @substep

        iris = RDatasets.dataset("datasets", "iris")
        scatter(Data(iris), Group(:Species), :SepalLength, :SepalWidth)
        @substep
        # use Position.stack to signal that you want bars stacked vertically rather than superimposed
        plot(Position.stack, histogram, Data(iris), Group(:Species), :SepalLength)
        @substep

        wireframe(
            density(trim=true),
            Data(iris), Group(:Species), :SepalLength, :SepalWidth,
            transparency = true, linewidth = 0.1
        )
        @substep

        scatter(
            Data(iris),
            Group(marker = :Species, color = bycolumn),
            :SepalLength, (:PetalLength, :PetalWidth)
        )
        @substep

        barplot(["hi", "ima", "string"], rand(3))
        @substep

    end
end


text: function div() end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
quoted(x) = headof(x) === :quote || headof(x) === :quotenode
unquoted(x) = isunarycall(x) && valof(x.args[1]) == "\$"

function remove_ref(x::EXPR)
    if hasref(x) && refof(x) isa Binding && refof(x).refs isa Vector
        for ia in enumerate(refof(x).refs)
            if ia[2] == x
                deleteat!(refof(x).refs, ia[1])
                setref!(x, nothing)
                return
            end
        end
        error()
    end
end

function clear_binding(x::EXPR)
    if bindingof(x) isa Binding
        for r in bindingof(x).refs
            if r isa EXPR
                setref!(r, nothing)
            elseif r isa Binding
                if r.type == bindingof(x)
                    r.type = nothing
                else
                    clear_binding(r)
                end
            end
        end
        x.meta.binding = nothing
    end
end
function clear_scope(x::EXPR)
    if hasmeta(x) && scopeof(x) isa Scope
        setparent!(scopeof(x), nothing)
        empty!(scopeof(x).names)
        if headof(x) === :file && scopeof(x).modules isa Dict && scopehasmodule(scopeof(x), :Base) && scopehasmodule(scopeof(x), :Core)
            m1, m2 = getscopemodule(scopeof(x), :Base), getscopemodule(scopeof(x), :Core)
            empty!(scopeof(x).modules)
            addmoduletoscope!(scopeof(x), m1)
            addmoduletoscope!(scopeof(x), m2)
        else
            scopeof(x).modules = nothing
        end
        if scopeof(x).overloaded !== nothing
            empty!(scopeof(x).overloaded)
        end
    end
end

function clear_ref(x::EXPR)
    if refof(x) isa Binding
        if refof(x).refs isa Vector
            for i in 1:length(refof(x).refs)
                if refof(x).refs[i] == x
                    deleteat!(refof(x).refs, i)
                    break
                end
            end
        end
        setref!(x, nothing)
    elseif refof(x) !== nothing
        setref!(x, nothing)
    end
end
function clear_error(x::EXPR)
    if hasmeta(x) && x.meta.error !== nothing
        x.meta.error = nothing
    end
end
function clear_meta(x::EXPR)
    clear_binding(x)
    clear_ref(x)
    clear_scope(x)
    clear_error(x)
    if x.args !== nothing
        for a in x.args
            clear_meta(a)
        end
    end
    # if x.trivia !== nothing
    #     for a in x.trivia
    #         clear_meta(a)
    #     end
    # end
end

function get_root_method(b, server)
    return b
end

function get_root_method(b::Binding, server)
    if CoreTypes.isfunction(b.type) && !isempty(b.refs)
        first(b.refs)
    else
        b
    end
end

function retrieve_delayed_scope(x)
    if (CSTParser.defines_function(x) || CSTParser.defines_macro(x)) && scopeof(x) !== nothing
        if parentof(scopeof(x)) !== nothing
            return parentof(scopeof(x))
        else
            return scopeof(x)
        end
    else
        return retrieve_scope(x)
    end
    return nothing
end

function retrieve_scope(x)
    if scopeof(x) !== nothing
        return scopeof(x)
    elseif parentof(x) isa EXPR
        return retrieve_scope(parentof(x))
    end
    return
end

# function find_return_statements(x::EXPR)
#     rets = EXPR[]
#     if CSTParser.defines_function(x)
#         find_return_statements(x.args[2], true, rets)
#     end
#     return rets
# end

# function find_return_statements(x::EXPR, last_stmt, rets)
#     if last_stmt && !(headof(x) === :block || headof(x) === :if || iskw(x))
#         push!(rets, x)
#         return rets, false
#     end

#     if headof(x) === :return
#         push!(rets, x)
#         return rets, true
#     end

#     for i = 1:length(x)
#         _, stop_iter = find_return_statements(x[i], last_stmt && (i == length(x) || (headof(x) === CSTParser.If && headof(x[i]) === CSTParser.Block)), rets)
#         stop_iter && break
#     end
#     return rets, false
# end

function find_exported_names(x::EXPR)
    exported_vars = EXPR[]
    for i in 1:length(x.args[3].args)
        expr = x.args[3].args[i]
        if headof(expr) === :export
            for j = 2:length(expr.args)
                if isidentifier(expr.args[j]) && hasref(expr.args[j])
                    push!(exported_vars, expr.args[j])
                end
            end
        end
    end
    return exported_vars
end

hasreadperm(p::String) = (uperm(p) & 0x04) == 0x04

# check whether a path is in (including subfolders) the julia base dir. Returns "" if not, and the path to the base dir if so.
function _is_in_basedir(path::String)
    i = findfirst(r".*base", path)
    i === nothing && return ""
    path1 = path[i]::String
    !hasreadperm(path1) && return ""
    !isdir(path1) && return ""
    files = readdir(path1)
    if all(f -> f in files, ["Base.jl", "coreio.jl", "essentials.jl", "exports.jl"])
        return path1
    end
    return ""
end

_is_macrocall_to_BaseDIR(arg) = headof(arg) === :macrocall && length(arg.args) == 2 && valof(arg.args[1]) == "@__DIR__"

isexportedby(k::Symbol, m::SymbolServer.ModuleStore) = haskey(m, k) && k in m.exportednames
isexportedby(k::String, m::SymbolServer.ModuleStore) = isexportedby(Symbol(k), m)
isexportedby(x::EXPR, m::SymbolServer.ModuleStore) = isexportedby(valof(x), m)
isexportedby(k, m::SymbolServer.ModuleStore) = false

function retrieve_toplevel_scope(x::EXPR)
    if scopeof(x) !== nothing && is_toplevel_scope(x)
        return scopeof(x)
    elseif parentof(x) isa EXPR
        return retrieve_toplevel_scope(parentof(x))
    else
        @info "Tried to reach toplevel scope, no scope found. Final expression $(headof(x))"
        return nothing
    end
end
retrieve_toplevel_scope(s::Scope) = (is_toplevel_scope(s) || !(parentof(s) isa Scope)) ? s : retrieve_toplevel_scope(parentof(s))
retrieve_toplevel_or_func_scope(s::Scope) = (is_toplevel_scope(s) || defines_function(s.expr) || !(parentof(s) isa Scope)) ? s : retrieve_toplevel_or_func_scope(parentof(s))

is_toplevel_scope(s::Scope) = is_toplevel_scope(s.expr)
is_toplevel_scope(x::EXPR) = CSTParser.defines_module(x) || headof(x) === :file

# b::SymbolServer.FunctionStore or DataTypeStore
# tls is a top-level Scope (expected to contain loaded modules)
# for a FunctionStore b, checks whether additional methods are provided by other packages
# f is a function that returns `true` if we want to break early from the loop

iterate_over_ss_methods(b, tls, env, f) = false
function iterate_over_ss_methods(b::SymbolServer.FunctionStore, tls::Scope, env::ExternalEnv, f)
    for m in b.methods
        ret = f(m)
        ret && return true
    end
    if b.extends in keys(getsymbolextendeds(env)) && tls.modules !== nothing
        # above should be modified,
        rootmod = SymbolServer._lookup(b.extends.parent, getsymbols(env)) # points to the module containing the initial function declaration
        if rootmod !== nothing && haskey(rootmod, b.extends.name) # check rootmod exists, and that it has the variable
            # find extensoions
            if haskey(getsymbolextendeds(env), b.extends) # method extensions listed
                for vr in getsymbolextendeds(env)[b.extends] # iterate over packages with extensions
                    !(SymbolServer.get_top_module(vr) in keys(tls.modules)) && continue
                    rootmod = SymbolServer._lookup(vr, getsymbols(env))
                    !(rootmod isa SymbolServer.ModuleStore) && continue
                    if haskey(rootmod.vals, b.extends.name) && (rootmod.vals[b.extends.name] isa SymbolServer.FunctionStore || rootmod.vals[b.extends.name] isa SymbolServer.DataTypeStore)# check package is available and has ref
                        for m in rootmod.vals[b.extends.name].methods #
                            ret = f(m)
                            ret && return true
                        end
                    end
                end
            end
        end
    end
    return false
end

function iterate_over_ss_methods(b::SymbolServer.DataTypeStore, tls::Scope, env::ExternalEnv, f)
    if b.name isa SymbolServer.VarRef
        bname = b.name
    elseif b.name isa SymbolServer.FakeTypeName
        bname = b.name.name
    end
    for m in b.methods
        ret = f(m)
        ret && return true
    end
    if (bname in keys(getsymbolextendeds(env))) && tls.modules !== nothing
        # above should be modified,
        rootmod = SymbolServer._lookup(bname.parent, getsymbols(env), true) # points to the module containing the initial function declaration
        if rootmod !== nothing && haskey(rootmod, bname.name) # check rootmod exists, and that it has the variable
            # find extensoions
            if haskey(getsymbolextendeds(env), bname) # method extensions listed
                for vr in getsymbolextendeds(env)[bname] # iterate over packages with extensions
                    !(SymbolServer.get_top_module(vr) in keys(tls.modules)) && continue
                    rootmod = SymbolServer._lookup(vr, getsymbols(env))
                    !(rootmod isa SymbolServer.ModuleStore) && continue
                    if haskey(rootmod.vals, bname.name) && (rootmod.vals[bname.name] isa SymbolServer.FunctionStore || rootmod.vals[bname.name] isa SymbolServer.DataTypeStore)# check package is available and has ref
                        for m in rootmod.vals[bname.name].methods #
                            ret = f(m)
                            ret && return true
                        end
                    end
                end
            end
        end
    end
    return false
end

"""
    is_in_fexpr(x::EXPR, f)
Check whether `x` isa the child of an expression for which `f(parent) == true`.
"""
is_in_fexpr(x::EXPR, f) = f(x) || (parentof(x) isa EXPR && is_in_fexpr(parentof(x), f))

"""
    get_in_fexpr(x::EXPR, f)
Get the `parent` of `x` for which `f(parent) == true`. (is_in_fexpr should be called first.)
"""
get_parent_fexpr(x::EXPR, f) = f(x) ? x : get_parent_fexpr(parentof(x), f)

maybe_get_parent_fexpr(x::Nothing, f) = nothing
maybe_get_parent_fexpr(x::EXPR, f) = f(x) ? x : maybe_get_parent_fexpr(parentof(x), f)

issigoffuncdecl(x::EXPR) = parentof(x) isa EXPR ? issigoffuncdecl(x, parentof(x)) : false
function issigoffuncdecl(x::EXPR, p::EXPR)
    if CSTParser.iswhere(p) || CSTParser.isdeclaration(p)
        return issigoffuncdecl(parentof(p))
    elseif CSTParser.defines_function(p)
        return true
    else
        return false
    end
end
issigoffuncdecl(x::EXPR, p) = false

function is_nameof_func(name)
    f = get_parent_fexpr(name, CSTParser.defines_function)
    f !== nothing && CSTParser.get_name(f) == name
end

function loose_refs(b::Binding)
    b.val isa EXPR || return b.refs # to account for `#global` binding which doesn't have a val
    scope = retrieve_scope(b.val)
    scope isa Scope && isidentifier(b.name) || return b.refs
    name_str = valofid(b.name)
    name_str isa String || return b.refs

    if is_soft_scope(scope) && parentof(scope) isa Scope && scopehasbinding(parentof(scope), name_str) && !scopehasbinding(scope, name_str)
        scope = parentof(scope)
    end
    state = LooseRefs(scope.expr, name_str, scope, [])
    state(scope.expr)
    vcat([r.refs for r in state.result]...)
end

mutable struct LooseRefs
    x::EXPR
    name::String
    scope::Scope
    result::Vector{Binding}
end

function (state::LooseRefs)(x::EXPR)
    if hasbinding(x)
        ex = bindingof(x).name
        if isidentifier(ex) && valofid(ex) == state.name
            push!(state.result, bindingof(x))
        end
    end
    if !hasscope(x) || (hasscope(x) && ((is_soft_scope(scopeof(x)) && !scopehasbinding(scopeof(x), state.name)) || scopeof(x) == state.scope))
        traverse(x, state)
    end
end


text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
using CartesianMesh
# using Revise
using HDGElasticity

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

function allequal(v1,v2)
    return all(v1 .== v2)
end

function allapprox(v1,v2,atol)
    @assert length(v1) == length(v2)
    return all([isapprox(v1[i],v2[i],atol=atol) for i = 1:length(v1)])
end

function plane_distance_function(coords,n,x0)
    return [n'*(coords[:,idx]-x0) for idx in 1:size(coords)[2]]
end

facequads = [[1,2],[],[1,2,3],[]]
solverid = HDGElasticity.face_to_solverid(facequads)
@test allequal(solverid,[1,0,2,0])

vbasis = TensorProductBasis(2,1)
mesh = UniformMesh([0.,0.],[1.,1.],[2,1])
coords = HDGElasticity.nodal_coordinates(mesh,vbasis)
NF = HDGElasticity.number_of_basis_functions(vbasis)
coeffs = reshape(plane_distance_function(coords,[1.,0.],[0.4,0.]),NF,:)
poly = InterpolatingPolynomial(1,vbasis)

dgmesh = HDGElasticity.DGMesh(mesh,coeffs,poly)
ufs = HDGElasticity.UniformFunctionSpace(dgmesh,1,4,coeffs,poly)

l1,m1 = (1.,2.)
l2,m2 = (2.,3.)
D1 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(l1,m1))
D2 = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(l2,m2))
stabilization = 1.0
cellmap = HDGElasticity.CellMap(dgmesh.domain[1])

solveridx = HDGElasticity.cell_to_solver_index(dgmesh.cellsign)
testsolveridx = [3 1
                 4 0]
@test allequal(solveridx,testsolveridx)

cellsolvers = HDGElasticity.CellSolvers(dgmesh,ufs,D1,D2,stabilization)
@test length(cellsolvers.localsolvers) == 4
@test allapprox(cellsolvers.stiffness[1],D1)
@test allapprox(cellsolvers.stiffness[2],D2)
@test cellsolvers.stabilization ≈ stabilization

update!(ufs.imap,ufs.icoeffs[1,1])
LL1 = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[1,1],ufs.fquads[1,1],
    dgmesh.facemaps,ufs.iquad,ufs.imap,-ufs.inormals[1],D1,stabilization,cellmap)
fLH1 = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,ufs.fquads[1,1],
    dgmesh.facemaps,ufs.fnormals,D1,stabilization,cellmap)
iLH1 = HDGElasticity.local_hybrid_operator_on_interface(ufs.vbasis,ufs.sbasis,
    ufs.iquad,ufs.imap,-ufs.inormals[1],D1,stabilization,cellmap)
LH1 = [hcat(fLH1...) iLH1]
@test allapprox(cellsolvers[1,1].LL,LL1)
@test allapprox(cellsolvers[1,1].LH,LH1)
@test allequal(cellsolvers[1,1].facetosolverid,[1,2,3,0,4])

LL2 = HDGElasticity.local_operator(ufs.vbasis,ufs.vquads[2,1],ufs.fquads[2,1],
    dgmesh.facemaps,ufs.iquad,ufs.imap,ufs.inormals[1],D2,stabilization,cellmap)
fLH2 = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,ufs.fquads[2,1],
    dgmesh.facemaps,ufs.fnormals,D2,stabilization,cellmap)
iLH2 = HDGElasticity.local_hybrid_operator_on_interface(ufs.vbasis,ufs.sbasis,
    ufs.iquad,ufs.imap,ufs.inormals[1],D2,stabilization,cellmap)
LH2 = [hcat(fLH2...) iLH2]
@test allapprox(cellsolvers[2,1].LL,LL2)
@test allapprox(cellsolvers[2,1].LH,LH2)
@test allapprox(cellsolvers[2,1].facetosolverid,[1,0,2,3,4])

uLL1 = HDGElasticity.local_operator(ufs.vbasis,ufs.vtpq,ufs.ftpq,
    dgmesh.facemaps,D1,stabilization,cellmap)
ufLH1 = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,ufs.ftpq,
    dgmesh.facemaps,ufs.fnormals,D1,stabilization,cellmap)
uLH1 = hcat(ufLH1...)
@test allapprox(cellsolvers.localsolvers[1].LL,uLL1)
@test allapprox(cellsolvers.localsolvers[1].LH,uLH1)
@test allequal(cellsolvers.localsolvers[1].facetosolverid,[1,2,3,4])

uLL2 = HDGElasticity.local_operator(ufs.vbasis,ufs.vtpq,ufs.ftpq,
    dgmesh.facemaps,D2,stabilization,cellmap)
ufLH2 = HDGElasticity.local_hybrid_operator(ufs.vbasis,ufs.sbasis,ufs.ftpq,
    dgmesh.facemaps,ufs.fnormals,D2,stabilization,cellmap)
uLH2 = hcat(ufLH2...)
@test allapprox(cellsolvers.localsolvers[2].LL,uLL2)
@test allapprox(cellsolvers.localsolvers[2].LH,uLH2)
@test allequal(cellsolvers.localsolvers[2].facetosolverid,[1,2,3,4])


text:
"""
	Constraint

	Abstract type to represent the different constraints of the problem

"""

abstract type Constraint end

"""
	isFiltrable

	Define if the contraint given in parameter is filtrable

"""

function isFiltrable end

"""
	filtrage

	Generic function that takes no parameter, modify the different bounds of the varaibles and return the modified variable

"""

function filtrage! end

include("Variables.jl")
include("Model.jl")
include("Change.jl")
#include("Constraints/AtMostOne.jl")
include("Constraints/Cardinalite.jl")
#include("Constraints/EmptyNIntersection.jl")
include("Constraints/EmptyIntersection.jl")
include("Constraints/Intersection.jl")
include("Constraints/DisjointUnion.jl")
include("Constraints/NDisjointUnion.jl")
include("Constraints/Social.jl")
include("Constraints/SortWeeks.jl")
include("Constraints/SortGroups.jl")


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ 8d375e12-5f3c-11ec-0551-0d5da3263f13
using DelimitedFiles

# ╔═╡ da0f707b-5c2a-4ab5-89c2-c48c04ef9404
d = readdlm("inputs/day_02.txt")

# ╔═╡ da7e5513-c0be-4bc6-9126-80ecfefd73db
directions = Dict(
	"forward" => [1, 0],
	"down" => [0, 1],
	"up" => [0, -1]
)

# ╔═╡ 127bd9b0-6b9b-4203-a325-3ddfae44dec2
instructions = map(command -> directions[command], d[:, 1])

# ╔═╡ 55369b8a-fd99-462f-a05a-6775a11652c3
steps = d[:, 2]

# ╔═╡ 4d5841f6-4f6f-44c3-af4a-084ce4e59864
*(sum(instructions .* steps)...)

# ╔═╡ 265dca33-2dd8-4adc-8612-4e905036ef5f
md"""
# Part II

In addition to horizontal position and depth, you'll also need to track a third value, aim, which also starts at 0. The commands also mean something entirely different than you first thought:

- down X increases your aim by X units.
- up X decreases your aim by X units.
- forward X does two things: It increases your horizontal position by X units. It increases your depth by your aim multiplied by X.

Again note that since you're on a submarine, down and up do the opposite of what you might expect: "down" means aiming in the positive direction.

"""

# ╔═╡ 8dca4413-94d8-4778-bda7-2143da76160b
function move_submarine(instructions, steps)
	aim = 0
	location = [0, 0]

	for (inst, step) in zip(instructions, steps)
		if inst == "up"
			aim -= step
		elseif inst == "down"
			aim += step
		else
			move = [step, step * aim]
			location += move
		end
	end
	return location
end

# ╔═╡ b9c2d28c-dc6c-49f9-9a6e-7ef2ba2dfec9
*(move_submarine(d[:, 1], steps)...)

# ╔═╡ 1b2e827f-31fc-4ba7-9d64-86ec46ce4ee8

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[[deps.DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[deps.Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"
"""

# ╔═╡ Cell order:
# ╠═8d375e12-5f3c-11ec-0551-0d5da3263f13
# ╠═da0f707b-5c2a-4ab5-89c2-c48c04ef9404
# ╠═da7e5513-c0be-4bc6-9126-80ecfefd73db
# ╠═127bd9b0-6b9b-4203-a325-3ddfae44dec2
# ╠═55369b8a-fd99-462f-a05a-6775a11652c3
# ╠═4d5841f6-4f6f-44c3-af4a-084ce4e59864
# ╠═265dca33-2dd8-4adc-8612-4e905036ef5f
# ╠═8dca4413-94d8-4778-bda7-2143da76160b
# ╠═b9c2d28c-dc6c-49f9-9a6e-7ef2ba2dfec9
# ╠═1b2e827f-31fc-4ba7-9d64-86ec46ce4ee8
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
module QMC

using ProgressMeter

using Measurements
using Statistics
using FFTW

using DelimitedFiles
using JLD2
using Printf

using DataStructures
using SparseArrays

import Base: zero
import Base: length, size, eltype, setindex!, getindex, firstindex, lastindex, rand, show

export BinaryQMCState, BinaryGroundState, BinaryThermalState,
        Hamiltonian, TFIM, nspins, nbonds, ClusterData, mc_step!, mc_step_beta!,
        resize_op_list!,
        sample, simulation_cell, magnetization, num_single_site_diag, num_single_site_offdiag,
        num_single_site, autocorrelation, correlation_time, jackknife, mean_and_stderr,
        lattice_bond_spins, ProbabilityHeap, ProbabilityVector, probability_vector

include("lattice.jl")
include("probabilityvector.jl")
include("op_sampler.jl")
include("hamiltonian.jl")
include("qmc_state.jl")
include("measurements.jl")
include("updates.jl")
include("error.jl")

end


text:
# Hestia is a registered Julia package
# Add Hestia.jl to your enviroment
# julia> ]
# pkg> add Hestia

using Hestia
λ = [10.0, 0.2, -0.0001]  # Thermal conductivity: λ(T) = 10 + 0.2 T - 1e-4 T^2
ρ = [7800.0]     # Mass density: constant
c = [330.0, 0.4] # Specific heat capacity: c(T) = 330 + 0.4 T

L = 0.3    # Length
W = 0.2    # Width
H = 0.1    # Height

Nx = 40     # Number of elements: x direction
Ny = 24     # Number of elements: y direction
Nz = 10
Ntotal = Nx*Ny*Nz

property = createDynamicIsoProperty(λ, ρ, c)
cuboid   = HeatCuboid(L, W, H,Nx, Ny, Nz, property)

### Boundaries ###
boundary = initBoundary(cuboid)

h = 10.0 # heat transfer coefficient
ϵ = 0.6  # emissivity
θamb = 300.0; # Ambient temperature
emission1 = createEmission(h, ϵ, θamb)  # Nonlinear BC: heat transfer (linear) and heat radiation (quartic/nonlinear)
emission2 = createEmission(h, 0.0, θamb)

setEmission!(boundary, emission1, :west) # Add emission to boundary
setEmission!(boundary, emission1, :east )
setEmission!(boundary, emission2, :north )
setEmission!(boundary, emission2, :topside)

# HeatProblem = container for geometry and boundary (not necessary)
# const heatproblem = CubicHeatProblem(cuboid, boundary)

## Simulation of cooling-down process ##
#=
cooling_down!(dθ, θ, param, t) = diffusion!(dθ, θ, cuboid, property, boundary)

# Initial conditions of ODE
θ₀ = 600.0  # Initial temperature
θinit = θ₀*ones(Ntotal)
tspan = (0.0, 200.0)
Δt    = 0.2             # Sampling time

import OrdinaryDiffEq
prob_cool = OrdinaryDiffEq.ODEProblem(cooling_down!,θinit,tspan)
sol_cool = OrdinaryDiffEq.solve(prob_cool,OrdinaryDiffEq.Euler(), dt=Δt, save_everystep=false)  # saveat=1.0)
=#

## Heating-up process ##
## Actuation ##
num_actuators = (4,3)        # Number of actuators
pos_actuators = :underside   # Position of actuators
actuation     = initActuation(cuboid)

# Create actuator characterization
scale     = 1.0;
power     = 3
curvature = 100.0;

config1  = setConfiguration(1.0, 3, 50)
config2  = setConfiguration(0.5, 2, 30)

setActuation!(actuation, cuboid, num_actuators, config1, :underside)

config_table = [config2 config1;
                config2 config1]
partition = [13 15;
             14  0]

#=
Partition on boundary east
x₃
^
| 15 |  0 |
|-------------
| 13 | 14 |
|------------> x₂
=#

setActuation!(actuation, cuboid, partition, config_table, :east)

num_act_total = num_actuators[1]*num_actuators[2] + 3
u_in = 4e5 * ones(num_act_total)

heating_up!(dθ, θ, param, t) = diffusion!(dθ, θ, cuboid, property, boundary, actuation, u_in)

## Simulation of heating-up process
# Initial conditions of ODE
θ₀ = 300.0  # Initial temperature
θinit = θ₀*ones(Ntotal)
tspan = (0.0, 200.0)
Δt    = 0.2  # Sampling time

import OrdinaryDiffEq

prob_heating = OrdinaryDiffEq.ODEProblem(heating_up!,θinit,tspan)
sol_heating  = OrdinaryDiffEq.solve(prob_heating,OrdinaryDiffEq.Euler(), dt=Δt, save_everystep=false)