text:
"""
    softmax(x; dims = 1)

[Softmax](https://en.wikipedia.org/wiki/Softmax_function) turns input array `x`
into probability distributions that sum to 1 along the dimensions specified by `dims`.
It is semantically equivalent to the following:

    softmax(x; dims = 1) = exp.(x) ./ sum(exp.(x), dims = dims)

with additional manipulations enhancing numerical stability.

For a matrix input `x` it will by default (`dims = 1`) treat it as a batch of vectors,
with each column independent. Keyword `dims = 2` will instead treat rows independently, and so on.

See also [`logsoftmax`](@ref).

# Examples

```jldoctest
julia> softmax([1, 2, 3])
3-element Vector{Float64}:
 0.09003057317038046
 0.24472847105479764
 0.6652409557748218

julia> softmax([1 2 3; 2 2 2])  # dims=1
2×3 Matrix{Float64}:
 0.268941  0.5  0.731059
 0.731059  0.5  0.268941

julia> softmax([1 2 3; 2 2 2]; dims=2)
2×3 Matrix{Float64}:
 0.0900306  0.244728  0.665241
 0.333333   0.333333  0.333333
```

Note that, when used with Flux.jl, `softmax` must not be passed to layers like `Dense`
which accept an activation function. The activation is broadcasted over the result,
thus applies to individual numbers. But `softmax` always needs to see the whole column.

```julia
julia> using Flux

julia> x = randn(Float32, 4, 4, 3, 13);

julia> model = Chain(Conv((4, 4), 3 => 8, tanh), Flux.flatten, Dense(8 => 7), softmax);

julia> model(x) |> size
(7, 13)

julia> Dense(4 => 7, softmax)(x)
ERROR: `softmax(x)` called with a number, but it expects an array.
```
"""
softmax(x::AbstractArray{T}; dims = 1) where {T} = softmax!(similar(x, float(T)), x; dims)

softmax!(x::AbstractArray; dims = 1) = softmax!(x, x; dims)

function softmax!(out::AbstractArray{T}, x::AbstractArray; dims = 1) where {T}
    max_ = maximum(x; dims)
    if all(isfinite, max_)
        @fastmath out .= exp.(x .- max_)
    else
        @fastmath @. out = ifelse(isequal(max_,Inf), ifelse(isequal(x,Inf), 1, 0), exp(x - max_))
    end
    out ./= sum(out; dims)
end

function ∇softmax_data(dy::AbstractArray{T}, y::AbstractArray{S}; dims = 1) where {T,S}
    dx = if within_grad()
        tmp = dy .* y
        tmp .- y .* sum(tmp; dims)
    else
        # This path is faster, only safe for 1st derivatives though.
        # Was previously `∇softmax!(dx, dy, x, y; dims)` to allow CUDA overloads,
        # but that was slow: https://github.com/FluxML/NNlibCUDA.jl/issues/30
        out = similar(y, promote_type(T,S))
        out .= dy .* y
        out .= out .- y .* sum(out; dims)
    end
end

function rrule(::typeof(softmax), x; dims = 1)
    y = softmax(x; dims)
    softmax_pullback(dy) = (NoTangent(), ∇softmax_data(unthunk(dy), y; dims))
    return y, softmax_pullback
end

within_grad() = false
rrule(::typeof(within_grad)) = true, _ -> (NoTangent(),)

"""
    logsoftmax(x; dims = 1)

Computes the log of softmax in a more numerically stable
way than directly taking `log.(softmax(xs))`. Commonly used in
computing cross entropy loss.

It is semantically equivalent to the following:

    logsoftmax(x; dims = 1) = x .- log.(sum(exp.(x), dims = dims))

See also [`softmax`](@ref).
"""
logsoftmax(x::AbstractArray{T}; dims = 1) where {T} = logsoftmax!(similar(x, float(T)), x; dims)

logsoftmax!(x::AbstractArray; dims = 1) = logsoftmax!(x, x; dims)

function logsoftmax!(out::AbstractArray{T}, x::AbstractArray; dims = 1) where {T}
    max_ = maximum(x; dims)
    if all(isfinite, max_)
        out .= x .- max_
    else
        @. out = ifelse(isequal(max_,Inf), ifelse(isequal(x,Inf), 0, -Inf), x - max_)
    end
    @fastmath log_ = log.(sum(exp, out; dims))
    out .-= log_
end

function ∇logsoftmax_data(dy::AbstractArray, y::AbstractArray; dims = 1)
    # This was previously `∇logsoftmax!(dx, dy, x, y; dims)` to allow CUDA overloads, but that was slow.
    dx = dy .- sum(dy; dims) .* exp.(y)
end

function rrule(::typeof(logsoftmax), x; dims = 1)
    y = logsoftmax(x; dims)
    logsoftmax_pullback(dy) = (NoTangent(), ∇logsoftmax_data(unthunk(dy), y; dims))
    return y, logsoftmax_pullback
end

"""
    logsumexp(x; dims = :)

Computes `log.(sum(exp.(x); dims))` in a numerically stable way.
Without `dims` keyword this returns a scalar.

See also [`logsoftmax`](@ref).
"""
function logsumexp(x::AbstractArray; dims = :)
    max_ = maximum(x; dims)
    @fastmath max_ .+ log.(sum(exp.(x .- max_); dims))
end

function rrule(::typeof(logsumexp), x; dims = :)
    # The gradient is `softmax`, but both compute `tmp` so it's worth saving.
    max_ = maximum(x; dims)
    @fastmath tmp = exp.(x .- max_)
    @fastmath y = max_ .+ log.(sum(tmp; dims))
    logsumexp_pullback(dy) = (NoTangent(), unthunk(dy) .* tmp ./ sum(tmp; dims))
    return y, logsumexp_pullback
end

# Informative error message if any of the softmax variants is called with a number
for f in (:softmax, :logsoftmax, :softmax!, :logsoftmax!, :logsumexp)
    @eval $(f)(x::Number, args...) =
      error("`", $(string(f)), "(x)` called with a number, but it expects an array. Usually this is because a layer like `Dense(3,4,softmax)` is broadcasting it like an activation function; `softmax` needs to be outside the layer.")
end


text:
using Plots

function cplot(f; interval = [-5,5], dx = 0.1, cmap = "hsv")
    x = interval[1]:dx:interval[2];
    X, Y = repmat(x',length(y),1), repmat(y,1,length(x));
    Z = f.(X + Y*1im);
    a = angle.(Z);
    pcolormesh(X, Y, a; cmap= cmap);
end

ext: jl

max_stars_repo_name: jeertmans/ELibraries

max_issues_repo_name: jeertmans/ELibraries

max_forks_repo_name: jeertmans/ELibraries


text:
# Mock

"""
    struct NewtonRaphsonSolver <:NonlinearSolver
      # Private fields
    end

Vanilla Newton-Raphson method
"""
struct NewtonRaphsonSolver <:NonlinearSolver
  ls::LinearSolver
  tol::Float64
  max_nliters::Int
end

struct NewtonRaphsonCache <: GridapType
  A::AbstractMatrix
  b::AbstractVector
  dx::AbstractVector
  ns::NumericalSetup
end

function solve!(x::AbstractVector,nls::NewtonRaphsonSolver,op::NonlinearOperator,cache::Nothing)
  b = residual(op, x)
  A = jacobian(op, x)
  dx = similar(b)
  ss = symbolic_setup(nls.ls, A)
  ns = numerical_setup(ss,A)
  _solve_nr!(x,A,b,dx,ns,nls,op)
  NewtonRaphsonCache(A,b,dx,ns)
end

function solve!(
  x::AbstractVector,nls::NewtonRaphsonSolver,op::NonlinearOperator,cache::NewtonRaphsonCache)
  b = cache.b
  A = cache.A
  dx = cache.dx
  ns = cache.ns
  residual!(b, op, x)
  jacobian!(A, op, x)
  numerical_setup!(ns,A)
  _solve_nr!(x,A,b,dx,ns,nls,op)
  cache
end

function _solve_nr!(x,A,b,dx,ns,nls,op)

  # Check for convergence on the initial residual
  isconv, conv0 = _check_convergence(nls,b)
  if isconv; return; end

  # Newton-like iterations
  for nliter in 1:nls.max_nliters

    # Solve linearized problem
    scale_entries!(b,-1)
    solve!(dx,ns,b)
    x .+= dx

    # Check convergence for the current residual
    residual!(b, op, x)
    isconv = _check_convergence(nls, b, conv0)
    if isconv; return; end

    if nliter == nls.max_nliters
      @unreachable
    end

    # Assemble jacobian (fast in-place version)
    # and prepare solver
    jacobian!(A, op, x)
    numerical_setup!(ns,A)

  end

end

function _check_convergence(nls,b)
  m0 = _inf_norm(b)
  (false, m0)
end

function _check_convergence(nls,b,m0)
  m = _inf_norm(b)
  m < nls.tol * m0
end

function _inf_norm(b)
  m = 0
  for bi in b
    m = max(m,abs(bi))
  end
  m
end

# Default concrete implementation

"""
    struct NLSolver <: NonlinearSolver
      # private fields
    end

The cache generated when using this solver has a field `result` that hosts the result
object generated by the underlying `nlsolve` function. It corresponds to the most latest solve.
"""
struct NLSolver <: NonlinearSolver
  ls::LinearSolver
  kwargs::Dict
end

"""
    NLSolver(ls::LinearSolver;kwargs...)
    NLSolver(;kwargs...)

Same kwargs as in `nlsolve`.
If `ls` is provided, it is not possible to use the `linsolve` kw-argument.
"""
function NLSolver(;kwargs...)
  ls = BackslashSolver()
  NLSolver(ls;kwargs...)
end

function NLSolver(ls::LinearSolver;kwargs...)
  @assert ! haskey(kwargs,:linsolve) "linsolve cannot be used here. It is managed internally"
  NLSolver(ls,kwargs)
end

mutable struct NLSolversCache <: GridapType
  f0::AbstractVector
  j0::AbstractMatrix
  df::OnceDifferentiable
  ns::NumericalSetup
  result
end

function solve!(x::AbstractVector,nls::NLSolver,op::NonlinearOperator,cache::Nothing)
  cache = _new_nlsolve_cache(x,nls,op)
  _nlsolve_with_updated_cache!(x,nls,op,cache)
  cache
end

function solve!(
  x::AbstractVector,nls::NLSolver,op::NonlinearOperator,cache::NLSolversCache)
  cache = _update_nlsolve_cache!(cache,x,op)
  _nlsolve_with_updated_cache!(x,nls,op,cache)
  cache
end

function _nlsolve_with_updated_cache!(x,nls,op,cache)
  df = cache.df
  ns = cache.ns
  kwargs = nls.kwargs
  function linsolve!(x,A,b)
    numerical_setup!(ns,A)
    solve!(x,ns,b)
  end
  r = nlsolve(df,x;linsolve=linsolve!,kwargs...)
  cache.result = r
  copy_entries!(x,r.zero)
end

function _new_nlsolve_cache(x0,nls,op)
  f!(r,x) = residual!(r,op,x)
  j!(j,x) = jacobian!(j,op,x)
  fj!(r,j,x) = residual_and_jacobian!(r,j,op,x)
  f0, j0 = residual_and_jacobian(op,x0)
  df = OnceDifferentiable(f!,j!,fj!,x0,f0,j0)
  ss = symbolic_setup(nls.ls,j0)
  ns = numerical_setup(ss,j0)
  NLSolversCache(f0,j0,df,ns,nothing)
end

function _update_nlsolve_cache!(cache,x0,op)
  f!(r,x) = residual!(r,op,x)
  j!(j,x) = jacobian!(j,op,x)
  fj!(r,j,x) = residual_and_jacobian!(r,j,op,x)
  f0 = cache.f0
  j0 = cache.j0
  ns = cache.ns
  residual_and_jacobian!(f0,j0,op,x0)
  df = OnceDifferentiable(f!,j!,fj!,x0,f0,j0)
  numerical_setup!(ns,j0)
  NLSolversCache(f0,j0,df,ns,nothing)
end


text:
##### Beginning of file

module Common # Begin submodule MirrorUpdater.Common

__precompile__(true)

import ..Types
import ..Utils

import ArgParse
import Dates
import HTTP
import Pkg
import TimeZones

function _toml_file_to_package(
        packagetoml_file_filename::String,
        )::Types.Package
    toml_file_text::String = read(packagetoml_file_filename, String)
    toml_file_parsed::Dict{String,Any}=Pkg.TOML.parse(toml_file_text)
    pkg_name::String = toml_file_parsed["name"]
    pkg_uuid::String = toml_file_parsed["uuid"]
    pkg_source_url::String = toml_file_parsed["repo"]
    pkg::Types.Package = Types.Package(
        ;
        name=pkg_name,
        uuid=pkg_uuid,
        source_url=pkg_source_url,
        )
    return pkg
end

function _get_uuid_from_toml_file(toml_file_filename::String)::String
    toml_file_text::String = read(toml_file_filename, String)
    toml_file_parsed::Dict{String,Any}=Pkg.TOML.parse(toml_file_text)
    uuid::String = toml_file_parsed["uuid"]
    return uuid
end

function _make_list(
        registry_list::Vector{Types.Registry},
        additional_repos::Vector{Types.SrcDestPair};
        do_not_try_url_list::Vector{String},
        try_but_allow_failures_url_list::Vector{String},
        )::Vector{Types.SrcDestPair}
    full_list::Vector{Types.SrcDestPair} = Types.SrcDestPair[]
    for x in additional_repos
        push!(full_list, x)
    end
    git = Utils._get_git_binary_path()
    for registry in registry_list
        registry_name = registry.name
        registry_uuid = registry.uuid
        registry_source_url = registry.url
        registry_destination_repo_name = _generate_destination_repo_name(
            registry
            )
        registry_src_dest_pair = Types.SrcDestPair(
            ;
            source_url = registry_source_url,
            destination_repo_name = registry_destination_repo_name,
            )
        push!(full_list, registry_src_dest_pair)
        if registry_source_url in do_not_try_url_list ||
                Types._name_with_git(registry_source_url) in do_not_try_url_list ||
                Types._name_without_git(registry_source_url) in do_not_try_url_list
            blah
            @warn(
                string(
                    "registry_source_url is in the do-not-try list, ",
                    "so skipping.",
                    ),
                registry_source_url,
                )
        else
            previous_dir::String = pwd()
            temp_dir_registry_git_clone_regular::String = mktempdir()
            cd(temp_dir_registry_git_clone_regular)
            cmd_git_clone_registry_regular = `$(git) clone $(registry.url)`
            @info(
                "Attempting to run command",
                cmd_git_clone_registry_regular,
                pwd(),
                ENV["PATH"],
                )
            clone_registry_regular_was_success =
                Utils.command_ran_successfully!!(
                    cmd_git_clone_registry_regular;
                    )
            if clone_registry_regular_was_success
                @info("Command ran successfully",)
                registry_toml_filename = joinpath(
                    temp_dir_registry_git_clone_regular,
                    registry_name,
                    "Registry.toml"
                    )
                registry_toml_file_uuid = _get_uuid_from_toml_file(
                    registry_toml_filename
                    )
                if lowercase(strip(registry_uuid)) !=
                        lowercase(strip(registry_toml_file_uuid))
                    error(
                        string(
                            "The UUID ($(registry_toml_file_uuid)) ",
                            "I found in the Registry.toml file does not ",
                            "match the UUID ($(registry_uuid)) ",
                            "that you provided.",
                            )
                        )
                end
                list_of_packagetoml_filenames::Vector{String} = String[]
                for (root, dirs, files) in
                        walkdir(temp_dir_registry_git_clone_regular)
                    for file in files
                        if lowercase(strip(file)) == "package.toml"
                            packagetoml_fn = joinpath(root, file)
                            push!(list_of_packagetoml_filenames, packagetoml_fn)
                        end
                    end
                end
                for packagetoml_file_filename in list_of_packagetoml_filenames
                    pkg = _toml_file_to_package(packagetoml_file_filename)
                    pkg_source_url = pkg.source_url
                    pkg_dest_repo_name = _generate_destination_repo_name(pkg)
                    pkg_src_dest_pair = Types.SrcDestPair(
                        ;
                        source_url=pkg_source_url,
                        destination_repo_name=pkg_dest_repo_name,
                        )
                    push!(full_list, pkg_src_dest_pair)
                end
            else
                @warn(
                    "Command did not run successfully",
                    cmd_git_clone_registry_regular,
                    pwd(),
                    ENV["PATH"],
                    )
                if registry_source_url in try_but_allow_failures_url_list ||
                        Types._name_with_git(registry_source_url) in
                            try_but_allow_failures_url_list ||
                        Types._name_without_git(registry_source_url) in
                            try_but_allow_failures_url_list
                    @warn(
                        string(
                            "URL is in the try-but-allow-failures list, ",
                            "so ignoring error ",
                            "that occured when running command",
                            ),
                        cmd_git_clone_registry_regular,
                        pwd(),
                        ENV["PATH"],
                        )
                else
                    error(
                        string(
                            "Encountered error when running command: ",
                            cmd_git_clone_registry_regular,
                            pwd(),
                            ENV["PATH"],
                            )
                        )
                end
            end
            cd(previous_dir)
            rm(
                temp_dir_registry_git_clone_regular;
                force = true,
                recursive = true,
                )
        end
    end
    unique_list_sorted::Vector{Types.SrcDestPair} = sort(unique(full_list))
    @info(
        string(
            "I made a list with ",
            "$(length(unique_list_sorted)) ",
            "unique pairs."
            )
        )
    return unique_list_sorted
end

function _generate_destination_repo_name(x::Types.Registry)::String
    result::String = string(
        strip(x.owner),
        "-",
        strip(x.name),
        "-",
        strip(x.uuid),
        )
    return result
end

function _generate_destination_repo_name(x::Types.Package)::String
    result::String = string(
        strip(x.name),
        "-",
        strip(x.uuid),
        )
    return result
end

function _src_dest_pair_to_string(x::Types.SrcDestPair)::String
    result::String = string(
        strip(x.source_url),
        " ",
        strip(x.destination_repo_name),
        )
    return result
end

function _src_dest_pair_list_to_string(
        v::Vector{Types.SrcDestPair}
        )::String
    v_sorted_unique::Vector{Types.SrcDestPair} = sort(unique(v))
    lines::Vector{String} = String[
        _src_dest_pair_to_string(x) for x in v_sorted_unique
        ]
    result::String = string(
        join(lines, "\n",),
        "\n",
        )
    return result
end

function _string_to_src_dest_pair_list(
        x::String
        )::Vector{Types.SrcDestPair}
    all_src_dest_pairs = Types.SrcDestPair[]
    lines::Vector{String} = convert(
        Vector{String},
        split(x, "\n",),
        )
    for line in lines
        columns::Vector{String} = convert(
            Vector{String},
            split(strip(line)),
            )
        if length(columns) == 2
            source_url::String = strip(columns[1])
            destination_repo_name::String = strip(columns[2])
            src_dest_pair::Types.SrcDestPair = Types.SrcDestPair(
                ;
                source_url = source_url,
                destination_repo_name = destination_repo_name,
                )
            push!(all_src_dest_pairs, src_dest_pair,)
        end
    end
    src_dest_pairs_sorted_unique::Vector{Types.SrcDestPair} = sort(
        unique(
            all_src_dest_pairs
            )
        )
    return src_dest_pairs_sorted_unique
end

function _remove_problematic_refs_before_github!!(
        ;
        packed_refs_filename::String,
        )::Nothing
    original_packed_refs_content::String = read(
        packed_refs_filename,
        String,
        )
    original_lines::Vector{String} = convert(
        Vector{String},
        split(strip(original_packed_refs_content), "\n")
        )
    function _line_is_ok_to_keep(x::String)::Bool
        result::Bool = (!(occursin("refs/pull/", x))) &&
            (!(occursin("cache/pull/", x)))
        return result
    end
    function _determine_new_gh_pages_branch_name(
            content::String,
            suggested_name::String,
            )::String
        if occursin(suggested_name, content)
            result = _determine_new_gh_pages_branch_name(
                content,
                string(suggested_name, "1"),
                )
        else
            result = suggested_name
        end
        return result
    end
    new_name_for_gh_pages_branch = _determine_new_gh_pages_branch_name(
        original_packed_refs_content,
        "gh-pages1",
        )
    function _transform_line(x::String)::String
        result_1::String = replace(
            x,
            "gh-pages" => new_name_for_gh_pages_branch,
            )
        return result_1
    end
    new_lines::Vector{String} = Vector{String}()
    for orig_line in original_lines
        if _line_is_ok_to_keep(orig_line)
            transformed_line = _transform_line(orig_line)
            push!(new_lines, transformed_line)
        end
    end
    new_packed_refs_content::String = string(
        join(new_lines, "\n"),
        "\n",
        )
    rm(
        packed_refs_filename;
        force = true,
        recursive = true,
        )
    write(packed_refs_filename, new_packed_refs_content)
    return nothing
end

function _push_mirrors!!(
        ;
        src_dest_pairs::Vector{Types.SrcDestPair},
        git_hosting_providers = git_hosting_providers,
        recursion_level::Integer = 0,
        max_recursion_depth::Integer = 10,
        is_dry_run::Bool = false,
        do_not_push_to_these_destinations::Vector{String},
        do_not_try_url_list::Vector{String},
        try_but_allow_failures_url_list::Vector{String},
        time_zone::Dates.TimeZone,
        )::Nothing
    @debug(string("Recursion level: $(recursion_level)"))
    git = Utils._get_git_binary_path()
    src_dest_pairs_sorted_unique::Vector{Types.SrcDestPair} = sort(
        unique(
            src_dest_pairs
            )
        )
    @info(
        string(
            "Running _push_mirrors!! with ",
            "$(length(src_dest_pairs_sorted_unique)) ",
            "unique pairs.",
            )
        )
    for pair_number = 1:length(src_dest_pairs_sorted_unique)
        @info(
            string(
                "Pair $(pair_number) of ",
                "$(length(src_dest_pairs_sorted_unique))",
                )
            )
        pair = src_dest_pairs_sorted_unique[pair_number]
        src_url = pair.source_url
        if src_url in do_not_try_url_list ||
                Types._name_with_git(src_url) in do_not_try_url_list ||
                Types._name_without_git(src_url) in do_not_try_url_list
            @warn(
                string("Src url is in the do not try list, so skipping."),
                src_url,
                )
        else
            previous_dir::String = pwd()
            temp_dir_repo_git_clone_regular::String = mktempdir()
            temp_dir_repo_git_clone_mirror::String = mktempdir()
            if recursion_level <= max_recursion_depth
                @info(
                    string(
                        "Now I will look for additional repos to mirror ",
                        "(e.g. BinaryBuilder repos that are referenced ",
                        "in this repo).",
                        )
                    )
                cd(temp_dir_repo_git_clone_regular)
                cmd_git_clone_repo_regular =
                    `$(git) clone $(src_url) GITCLONEREPOREGULAR`
                @info(
                    "Attempting to run command",
                    cmd_git_clone_repo_regular,
                    pwd(),
                    ENV["PATH"],
                    )
                repo_regular_clone_was_success =
                    Utils.command_ran_successfully!!(
                    cmd_git_clone_repo_regular;
                    )
                if repo_regular_clone_was_success
                    @info("Command ran successfully",)
                    cd(
                        joinpath(
                            temp_dir_repo_git_clone_regular,
                            "GITCLONEREPOREGULAR",
                            )
                        )
                    git_grep_results::String = try
                        strip(read(`$(git) grep Builder`, String))
                    catch exception
                        @info("ignoring exception: ", exception)
                        ""
                    end
                    list_of_new_src_dest_pairs::Vector{Types.SrcDestPair} =
                        Types.SrcDestPair[]
                    if length(git_grep_results) > 0
                        bin_bldr_pair_list::Vector{Types.SrcDestPair} =
                            _get_list_of_binary_builder_repos(
                                git_grep_results
                                )
                        for bin_bldr_pair in bin_bldr_pair_list
                            if bin_bldr_pair in src_dest_pairs
                            else
                                push!(
                                    list_of_new_src_dest_pairs,
                                    bin_bldr_pair,
                                    )
                            end
                        end
                    end
                    if (length(git_grep_results) > 0) &&
                            (length(list_of_new_src_dest_pairs) > 0)
                        if length(list_of_new_src_dest_pairs) == 1
                            @info(
                                string(
                                    "I found ",
                                    "1 ",
                                    "additional repo to mirror. ",
                                    "I will mirror ",
                                    " it first, and then I will return ",
                                    "to my previous list.",
                                    )
                                )
                        else
                            @info(
                                string(
                                    "I found ",
                                    "$(length(list_of_new_src_dest_pairs)) ",
                                    "additional repos to mirror. ",
                                    "I will mirror ",
                                    " them first, and then I will return ",
                                    "to my previous list.",
                                    )
                                )
                        end
                        _push_mirrors!!(
                            ;
                            do_not_push_to_these_destinations =
                                do_not_push_to_these_destinations,
                            src_dest_pairs = list_of_new_src_dest_pairs,
                            git_hosting_providers = git_hosting_providers,
                            recursion_level = recursion_level + 1,
                            max_recursion_depth = max_recursion_depth,
                            is_dry_run = is_dry_run,
                            do_not_try_url_list = do_not_try_url_list,
                            time_zone = time_zone,
                            try_but_allow_failures_url_list =
                                try_but_allow_failures_url_list,
                            )
                    else
                        @info(
                            string(
                                "I did not find any additional ",
                                "repos to mirror.",
                                )
                            )
                    end
                else
                    if src_url in try_but_allow_failures_url_list ||
                            Types._name_with_git(src_url) in
                                try_but_allow_failures_url_list ||
                            Types._name_without_git(src_url) in
                                try_but_allow_failures_url_list
                        @warn(
                            string(
                                "URL in the try-but-allow-failures list, ",
                                "so ignoring the error ",
                                "that occured while running command",
                                ),
                            cmd_git_clone_repo_regular,
                            pwd(),
                            ENV["PATH"],
                            )
                    else
                        error(
                            string(
                                "Encountered error when running command: ",
                                cmd_git_clone_repo_regular,
                                pwd(),
                                ENV["PATH"],
                                )
                            )
                    end
                end
            else
                @warn(
                    string(
                        "I have exceeded the maximum recursion depth.",
                        ),
                    recursion_level,
                    max_recursion_depth,
                    )
            end
            cd(temp_dir_repo_git_clone_mirror)
            cmd_git_repo_clone_mirror =
                `$(git) clone --mirror $(src_url) GITCLONEREPOMIRROR`
            @info(
                "Attempting to run command",
                cmd_git_repo_clone_mirror,
                pwd(),
                ENV["PATH"],
                )
            repo_mirror_clone_was_success =
                Utils.command_ran_successfully!!(
                cmd_git_repo_clone_mirror;
                )
            if repo_mirror_clone_was_success
                @info("Command ran successfully",)
                cd(
                    joinpath(
                        temp_dir_repo_git_clone_mirror,
                        "GITCLONEREPOMIRROR",
                        )
                    )
                @info("Processing the repository")
                packed_refs_filename = joinpath(
                    temp_dir_repo_git_clone_mirror,
                    "GITCLONEREPOMIRROR",
                    "packed-refs",
                    )
                _remove_problematic_refs_before_github!!(
                    ;
                    packed_refs_filename = packed_refs_filename,
                    )
                destination_repo_name = pair.destination_repo_name
                if destination_repo_name in do_not_push_to_these_destinations ||
                        destination_repo_name in do_not_push_to_these_destinations ||
                        destination_repo_name in do_not_push_to_these_destinations
                    @warn(
                        string(
                            "Destination repo name is in the ",
                            "do_not_push_to_these_destinations list, ",
                            "so skipping.",
                            )
                        )
                else
                    if is_dry_run
                        @info(
                            string(
                                "This is a dry run, so I will not ",
                                "push to any git hosting providers.",
                                )
                            )
                    else
                        for p = 1:length(git_hosting_providers)
                            @info(
                                string(
                                    "Git hosting provider ",
                                    "$(p) of ",
                                    "$(length(git_hosting_providers))",
                                    ),
                                )
                            provider = git_hosting_providers[p]
                            args1 = Dict(
                                :repo_name => destination_repo_name,
                                )
                            @info(
                                string(
                                    "Making sure that repo exists on ",
                                    "git hosting provider $(p). ",
                                    "(If it does not already exist, ",
                                    "I will create it.)",
                                    )
                                )
                            provider(:create_repo)(args1)
                            args2 = Dict(
                                :repo_name => destination_repo_name,
                                :directory => pwd(),
                                :git => git,
                                :try_but_allow_failures_url_list =>
                                    try_but_allow_failures_url_list,
                                )
                            @info(
                                string(
                                    "Attempting to push to ",
                                    "git hosting provider $(p).",
                                    )
                                )
                            push_to_provider_was_success = try
                                provider(:push_mirrored_repo)(args2)
                                true
                            catch exception
                                @warn(
                                    "ignoring exception: ",
                                    exception,
                                    )
                                false
                            end
                            if push_to_provider_was_success
                                when = Dates.now(TimeZones.localzone(),)
                                args3 = Dict(
                                    :source_url => src_url,
                                    :when => when,
                                    :time_zone => time_zone,
                                    )
                                repo_description_default::String =
                                    Utils.default_repo_description(
                                        ;
                                        from = src_url,
                                        when = when,
                                        time_zone = time_zone,
                                        )
                                repo_description_provider::String = try
                                    provider(
                                        :generate_new_repo_description)(
                                        args3)
                                catch exception
                                    @warn(
                                        string("ignoring exception: "),
                                        exception,
                                        )
                                    ""
                                end

                                new_repo_description::String = ""
                                if length(
                                        strip(
                                            repo_description_provider
                                            )
                                        ) == 0
                                    new_repo_description = strip(
                                        repo_description_default
                                        )
                                else
                                    new_repo_description = strip(
                                        repo_description_provider
                                        )
                                end

                                @debug(
                                    string("Repo descriptions: "),
                                    repo_description_default,
                                    repo_description_provider,
                                    new_repo_description,
                                    )

                                args4 = Dict(
                                    :repo_name =>
                                        destination_repo_name,
                                    :new_repo_description =>
                                        new_repo_description
                                    )
                                @info(
                                    string(
                                        "Attempting to update ",
                                        "repo description on git hosting ",
                                        "provider $(p).",
                                        ),
                                    destination_repo_name,
                                    new_repo_description,
                                    )
                                provider(:update_repo_description)(args4)
                            else
                                error(
                                    string(
                                        "Push to provider $(p) ",
                                        "was not a success.",
                                        )
                                    )
                            end
                        end
                    end
                end
            else
                if src_url in try_but_allow_failures_url_list ||
                        Types._name_with_git(src_url) in
                            try_but_allow_failures_url_list ||
                        Types._name_without_git(src_url) in
                            try_but_allow_failures_url_list
                    @warn(
                        string(
                            "URL in the try-but-allow-failures list, ",
                            "so ignoring the error ",
                            "that occured while running command",
                            ),
                        cmd_git_repo_clone_mirror,
                        pwd(),
                        ENV["PATH"],
                        )
                else
                    error(
                        string(
                            "Encountered error when running command: ",
                            cmd_git_repo_clone_mirror,
                            pwd(),
                            ENV["PATH"],
                            )
                        )
                end
            end
            cd(previous_dir)
            rm(
                temp_dir_repo_git_clone_regular;
                force = true,
                recursive = true,
                )
            rm(
                temp_dir_repo_git_clone_mirror;
                force = true,
                recursive = true,
                )
        end
    end
    return nothing
end

function _get_list_of_binary_builder_repos(
        text::AbstractString,
        )::Vector{Types.SrcDestPair}
    result::Vector{Types.SrcDestPair} = Types.SrcDestPair[]
    lines::Vector{String} = convert(
        Vector{String},
        split(strip(text), "\n"),
        )
    regex_1::Regex = r"https:\/\/github.com\/(\w*?)\/(\w*?)\/"
    for line in lines
        line_stripped::String = strip(line)
        if occursin(regex_1, line_stripped)
            regex_match::RegexMatch = match(regex_1, line_stripped)
            github_repo_owner::String = strip(regex_match[1])
            github_repo_name::String = strip(regex_match[2])
            source_url::String = string(
                "https://github.com/",
                github_repo_owner,
                "/",
                github_repo_name,
                )
            destination_repo_name::String = string(
                github_repo_owner,
                "-",
                github_repo_name,
                )
            new_pair = Types.SrcDestPair(
                ;
                source_url = source_url,
                destination_repo_name = destination_repo_name,
                )
            push!(
                result,
                new_pair,
                )
        else
        end
    end
    return result
end

function _add_trailing_spaces(x::AbstractString, n::Integer)::String
    temp::String = strip(convert(String, x))
    if length(temp) >= n
        result::String = temp
    else
        deficit::Int = n - length(temp)
        result = string(temp, repeat(" ", deficit))
    end
    return result
end

function _interval_contains_x(
        interval::Types.AbstractInterval,
        pair::Types.SrcDestPair,
        )::Bool
    result::Bool = _interval_contains_x(
        interval,
        pair.destination_repo_name,
        )
    return result
end

function _pairs_that_fall_in_interval(
        list_of_pairs::Vector{Types.SrcDestPair},
        interval::Types.AbstractInterval,
        )::Vector{Types.SrcDestPair}
    ith_pair_falls_in_interval::Vector{Bool} = Vector{Bool}(
        undef,
        length(list_of_pairs),
        )
    for i = 1:length(list_of_pairs)
        ith_pair = list_of_pairs[i]
        ith_pair_falls_in_interval[i] = _interval_contains_x(
            interval,
            ith_pair,
            )
    end
    full_sublist::Vector{Types.SrcDestPair} = list_of_pairs[
        ith_pair_falls_in_interval
        ]
    unique_sorted_sublist::Vector{Types.SrcDestPair} = sort(
        unique(full_sublist)
        )
    return unique_sorted_sublist
end

function _interval_contains_x(
        interval::Types.NoBoundsInterval,
        x::AbstractString,
        )::Bool
    result::Bool = true
    return result
end

function _interval_contains_x(
        interval::Types.LowerAndUpperBoundInterval,
        x::AbstractString,
        )::Bool
    x_stripped::String = strip(convert(String, x))
    left::String = strip(interval.left)
    right::String = strip(interval.right)
    result::Bool = (left <= x_stripped) && (x_stripped < right)
    return result
end

function _interval_contains_x(
        interval::Types.LowerBoundOnlyInterval,
        x::AbstractString,
        )::Bool
    x_stripped::String = strip(convert(String, x))
    left::String = strip(interval.left)
    result::Bool = left <= x_stripped
    return result
end

function _interval_contains_x(
        interval::Types.UpperBoundOnlyInterval,
        x::AbstractString,
        )
    x_stripped::String = strip(convert(String, x))
    right::String = strip(interval.right)
    result::Bool = x_stripped < right
    return result
end

end # End submodule MirrorUpdater.Common

##### End of file


text:
using SIMeasurements: Meter, Kilogram, unit_symbol, dimensionality
import SIMeasurements: conversion_factors
facts("Basic SI Types") do
    @fact dimensionality(Meter) --> [1, 0, 0, 0, 0, 0, 0, 0, 0]
    @fact unit_symbol(Meter) --> "m"

    @fact dimensionality(Kilogram) --> [0, 1, 0, 0, 0, 0, 0, 0, 0]
    @fact unit_symbol(Kilogram) --> "kg"
end

facts("Unit symbols of more complicated types") do
    const Kludge = Unit{:SI, Dimensions(luminosity=-1, mass=2, length=3)}()
    @fact dimensionality(Kludge) --> [3, 2, 0, 0, 0, 0, -1, 0, 0]
    @fact unit_symbol(Kludge) --> "m³⋅kg²⋅Cd⁻¹"
end

facts("Operations over units") do
    @fact dimensionality(Meter * Meter) --> [2, 0, 0, 0, 0, 0, 0, 0, 0]
    @fact typeof(dimensionality(Meter * Meter)) --> exactly(Dimensions)
    @fact dimensionality(Meter * Kilogram) --> [1, 1, 0, 0, 0, 0, 0, 0, 0]
    @fact typeof(dimensionality(Meter * Kilogram)) --> exactly(Dimensions)
    @fact dimensionality(Meter / Kilogram) --> [1, -1, 0, 0, 0, 0, 0, 0, 0]
    @fact typeof(dimensionality(Meter / Kilogram)) --> exactly(Dimensions)
    @fact dimensionality(Meter / Meter) --> [0, 0, 0, 0, 0, 0, 0, 0, 0]
    @fact typeof(dimensionality(Meter / Meter)) --> exactly(Dimensions)

    @fact dimensionality(Meter^-1) --> [-1, 0, 0, 0, 0, 0, 0, 0, 0]
    @fact typeof(dimensionality(Meter^-1)) --> exactly(Dimensions)
    @fact dimensionality(Meter^0) --> [0, 0, 0, 0, 0, 0, 0, 0, 0]
    @fact typeof(dimensionality(Meter^0)) --> exactly(Dimensions)
end

const Kloodge = Unit{:Beware, Dimensions(length=1)}()
conversion_factors(::typeof(Kloodge), ::typeof(Meter)) = 3//10, 3//5
const Blagger = Unit{:Beware, Dimensions(mass=2)}()
conversion_factors(::Unit{:Beware, Dimensions(mass=1)}, ::typeof(Kilogram)) =
    2//3, 0
facts("Conversion factor") do
   @fact conversion_factors(Kloodge, Meter) --> (3//10, 3//5)
   @fact conversion_factors(Meter, Kloodge) --> (10//3, -2)
   @fact conversion_factors(Meter^2, Kloodge^2) --> (100//9, 0)
   @fact conversion_factors(Blagger, Kilogram^2) --> (4//9, 0)
end


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

# TODO: Taal refactor
# - analysis_interval part as PeriodicCallback called after a certain amount of simulation time
"""
    AnalysisCallback(semi; interval=0,
                           save_analysis=false,
                           output_directory="out",
                           analysis_filename="analysis.dat",
                           extra_analysis_errors=Symbol[],
                           extra_analysis_integrals=())

Analyze a numerical solution every `interval` time steps and print the
results to the screen. If `save_analysis`, the results are also saved in
`joinpath(output_directory, analysis_filename)`.

Additional errors can be computed, e.g. by passing `extra_analysis_errors = [:primitive]`.

Further scalar functions `func` in `extra_analysis_integrals` are applied to the numerical
solution and integrated over the computational domain.
See `Trixi.analyze`, `Trixi.pretty_form_utf`, `Trixi.pretty_form_ascii` for further
information on how to create custom analysis quantities.
"""
mutable struct AnalysisCallback{Analyzer, AnalysisIntegrals, InitialStateIntegrals, Cache}
  start_time::Float64
  interval::Int
  save_analysis::Bool
  output_directory::String
  analysis_filename::String
  analyzer::Analyzer
  analysis_errors::Vector{Symbol}
  analysis_integrals::AnalysisIntegrals
  initial_state_integrals::InitialStateIntegrals
  cache::Cache
end

# TODO: Taal bikeshedding, implement a method with less information and the signature
# function Base.show(io::IO, analysis_callback::AnalysisCallback)
# end
function Base.show(io::IO, ::MIME"text/plain", cb::DiscreteCallback{<:Any, <:AnalysisCallback})
  @nospecialize cb # reduce precompilation time

  if get(io, :compact, false)
    show(io, cb)
  else
    analysis_callback = cb.affect!

    setup = Pair{String,Any}[
             "interval" => analysis_callback.interval,
             "analyzer" => analysis_callback.analyzer,
            ]
    for (idx, error) in enumerate(analysis_callback.analysis_errors)
      push!(setup, "│ error " * string(idx) => error)
    end
    for (idx, integral) in enumerate(analysis_callback.analysis_integrals)
      push!(setup, "│ integral " * string(idx) => integral)
    end
    push!(setup, "save analysis to file" => analysis_callback.save_analysis ? "yes" : "no")
    if analysis_callback.save_analysis
      push!(setup, "│ filename" => analysis_callback.analysis_filename)
      push!(setup, "│ output directory" => abspath(normpath(analysis_callback.output_directory)))
    end
    summary_box(io, "AnalysisCallback", setup)
  end
end

function AnalysisCallback(semi::AbstractSemidiscretization; kwargs...)
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  AnalysisCallback(mesh, equations, solver, cache; kwargs...)
end

function AnalysisCallback(mesh, equations::AbstractEquations, solver, cache;
                          interval=0,
                          save_analysis=false,
                          output_directory="out",
                          analysis_filename="analysis.dat",
                          extra_analysis_errors=Symbol[],
                          analysis_errors=union(default_analysis_errors(equations), extra_analysis_errors),
                          extra_analysis_integrals=(),
                          analysis_integrals=union(default_analysis_integrals(equations), extra_analysis_integrals),
                          RealT=real(solver),
                          uEltype=eltype(cache.elements),
                          kwargs...)
  # Decide when the callback is activated.
  # With error-based step size control, some steps can be rejected. Thus,
  #   `integrator.iter >= integrator.destats.naccept`
  #    (total #steps)       (#accepted steps)
  # We need to check the number of accepted steps since callbacks are not
  # activated after a rejected step.
  condition = (u, t, integrator) -> interval > 0 && ( (integrator.destats.naccept % interval == 0 &&
                                                       !(integrator.destats.naccept == 0 && integrator.iter > 0)) ||
                                                     isfinished(integrator))

  analyzer = SolutionAnalyzer(solver; kwargs...)
  cache_analysis = create_cache_analysis(analyzer, mesh, equations, solver, cache, RealT, uEltype)

  analysis_callback = AnalysisCallback(0.0, interval, save_analysis, output_directory, analysis_filename,
                                       analyzer,
                                       analysis_errors, Tuple(analysis_integrals),
                                       SVector(ntuple(_ -> zero(uEltype), Val(nvariables(equations)))),
                                       cache_analysis)

  DiscreteCallback(condition, analysis_callback,
                   save_positions=(false,false),
                   initialize=initialize!)
end

function initialize!(cb::DiscreteCallback{Condition,Affect!}, u_ode, t, integrator) where {Condition, Affect!<:AnalysisCallback}
  semi = integrator.p
  initial_state_integrals = integrate(u_ode, semi)
  _, equations, _, _ = mesh_equations_solver_cache(semi)

  analysis_callback = cb.affect!
  analysis_callback.initial_state_integrals = initial_state_integrals
  @unpack save_analysis, output_directory, analysis_filename, analysis_errors, analysis_integrals = analysis_callback

  if save_analysis && mpi_isroot()
    mkpath(output_directory)

    # write header of output file
    open(joinpath(output_directory, analysis_filename), "w") do io
      @printf(io, "#%-8s", "timestep")
      @printf(io, "  %-14s", "time")
      @printf(io, "  %-14s", "dt")
      if :l2_error in analysis_errors
        for v in varnames(cons2cons, equations)
          @printf(io, "   %-14s", "l2_" * v)
        end
      end
      if :linf_error in analysis_errors
        for v in varnames(cons2cons, equations)
          @printf(io, "   %-14s", "linf_" * v)
        end
      end
      if :conservation_error in analysis_errors
        for v in varnames(cons2cons, equations)
          @printf(io, "   %-14s", "cons_" * v)
        end
      end
      if :residual in analysis_errors
        for v in varnames(cons2cons, equations)
          @printf(io, "   %-14s", "res_" * v)
        end
      end
      if :l2_error_primitive in analysis_errors
        for v in varnames(cons2prim, equations)
          @printf(io, "   %-14s", "l2_" * v)
        end
      end
      if :linf_error_primitive in analysis_errors
        for v in varnames(cons2prim, equations)
          @printf(io, "   %-14s", "linf_" * v)
        end
      end

      for quantity in analysis_integrals
        @printf(io, "   %-14s", pretty_form_ascii(quantity))
      end

      println(io)
    end

  end

  analysis_callback.start_time = time_ns()
  analysis_callback(integrator)
  return nothing
end

# TODO: Taal refactor, allow passing an IO object (which could be devnull to avoid cluttering the console)
function (analysis_callback::AnalysisCallback)(integrator)
  semi = integrator.p
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  @unpack dt, t = integrator
  iter = integrator.destats.naccept

  runtime_absolute = 1.0e-9 * (time_ns() - analysis_callback.start_time)
  runtime_relative = 1.0e-9 * take!(semi.performance_counter) / ndofs(semi)

  @trixi_timeit timer() "analyze solution" begin
    # General information
    mpi_println()
    mpi_println("─"^100)
    # TODO: Taal refactor, polydeg is specific to DGSEM
    mpi_println(" Simulation running '", get_name(equations), "' with ", summary(solver))
    mpi_println("─"^100)
    mpi_println(" #timesteps:     " * @sprintf("% 14d", iter) *
                "               " *
                " run time:       " * @sprintf("%10.8e s", runtime_absolute))
    mpi_println(" Δt:             " * @sprintf("%10.8e", dt) *
                "               " *
                " time/DOF/rhs!:  " * @sprintf("%10.8e s", runtime_relative))
    mpi_println(" sim. time:      " * @sprintf("%10.8e", t))
    mpi_println(" #DOF:           " * @sprintf("% 14d", ndofs(semi)))
    mpi_println(" #elements:      " * @sprintf("% 14d", nelements(mesh, solver, cache)))

    # Level information (only show for AMR)
    print_amr_information(integrator.opts.callback, mesh, solver, cache)
    mpi_println()

    # Open file for appending and store time step and time information
    if mpi_isroot() && analysis_callback.save_analysis
      io = open(joinpath(analysis_callback.output_directory, analysis_callback.analysis_filename), "a")
      @printf(io, "% 9d", iter)
      @printf(io, "  %10.8e", t)
      @printf(io, "  %10.8e", dt)
    else
      io = devnull
    end

    # Calculate current time derivative (needed for semidiscrete entropy time derivative, residual, etc.)
    du_ode = first(get_tmp_cache(integrator))
    @notimeit timer() rhs!(du_ode, integrator.u, semi, t)
    u  = wrap_array(integrator.u, mesh, equations, solver, cache)
    du = wrap_array(du_ode,       mesh, equations, solver, cache)
    l2_error, linf_error = analysis_callback(io, du, u, integrator.u, t, semi)

    mpi_println("─"^100)
    mpi_println()

    # Add line break and close analysis file if it was opened
    if mpi_isroot() && analysis_callback.save_analysis
      # This resolves a possible type instability introduced above, since `io`
      # can either be an `IOStream` or `devnull`, but we know that it must be
      # an `IOStream here`.
      println(io::IOStream)
      close(io::IOStream)
    end
  end

  # avoid re-evaluating possible FSAL stages
  u_modified!(integrator, false)

  # Return errors for EOC analysis
  return l2_error, linf_error
end

# This method is just called internally from `(analysis_callback::AnalysisCallback)(integrator)`
# and serves as a function barrier. Additionally, it makes the code easier to profile and optimize.
function (analysis_callback::AnalysisCallback)(io, du, u, u_ode, t, semi)
  @unpack analyzer, analysis_errors, analysis_integrals = analysis_callback
  cache_analysis = analysis_callback.cache
  _, equations, _, _ = mesh_equations_solver_cache(semi)

  # Calculate and print derived quantities (error norms, entropy etc.)
  # Variable names required for L2 error, Linf error, and conservation error
  if any(q in analysis_errors for q in
         (:l2_error, :linf_error, :conservation_error, :residual)) && mpi_isroot()
    print(" Variable:    ")
    for v in eachvariable(equations)
      @printf("   %-14s", varnames(cons2cons, equations)[v])
    end
    println()
  end

  # Calculate L2/Linf errors, which are also returned
  l2_error, linf_error = calc_error_norms(u_ode, t, analyzer, semi, cache_analysis)

  if mpi_isroot()
    # L2 error
    if :l2_error in analysis_errors
      print(" L2 error:    ")
      for v in eachvariable(equations)
        @printf("  % 10.8e", l2_error[v])
        @printf(io, "  % 10.8e", l2_error[v])
      end
      println()
    end

    # Linf error
    if :linf_error in analysis_errors
      print(" Linf error:  ")
      for v in eachvariable(equations)
        @printf("  % 10.8e", linf_error[v])
        @printf(io, "  % 10.8e", linf_error[v])
      end
      println()
    end
  end

  # Conservation errror
  if :conservation_error in analysis_errors
    @unpack initial_state_integrals = analysis_callback
    state_integrals = integrate(u_ode, semi)

    if mpi_isroot()
      print(" |∑U - ∑U₀|:  ")
      for v in eachvariable(equations)
        err = abs(state_integrals[v] - initial_state_integrals[v])
        @printf("  % 10.8e", err)
        @printf(io, "  % 10.8e", err)
      end
      println()
    end
  end

  # Residual (defined here as the vector maximum of the absolute values of the time derivatives)
  if :residual in analysis_errors
    mpi_print(" max(|Uₜ|):   ")
    for v in eachvariable(equations)
      # Calculate maximum absolute value of Uₜ
      res = maximum(abs, view(du, v, ..))
      if mpi_isparallel()
        # TODO: Debugging, here is a type instability
        global_res = MPI.Reduce!(Ref(res), max, mpi_root(), mpi_comm())
        if mpi_isroot()
          res::eltype(du) = global_res[]
        end
      end
      if mpi_isroot()
        @printf("  % 10.8e", res)
        @printf(io, "  % 10.8e", res)
      end
    end
    mpi_println()
  end

  # L2/L∞ errors of the primitive variables
  if :l2_error_primitive in analysis_errors || :linf_error_primitive in analysis_errors
    l2_error_prim, linf_error_prim = calc_error_norms(cons2prim, u_ode, t, analyzer, semi, cache_analysis)

    if mpi_isroot()
      print(" Variable:    ")
      for v in eachvariable(equations)
        @printf("   %-14s", varnames(cons2prim, equations)[v])
      end
      println()

      # L2 error
      if :l2_error_primitive in analysis_errors
        print(" L2 error prim.: ")
        for v in eachvariable(equations)
          @printf("%10.8e   ", l2_error_prim[v])
          @printf(io, "  % 10.8e", l2_error_prim[v])
        end
        println()
      end

      # L∞ error
      if :linf_error_primitive in analysis_errors
        print(" Linf error pri.:")
        for v in eachvariable(equations)
          @printf("%10.8e   ", linf_error_prim[v])
          @printf(io, "  % 10.8e", linf_error_prim[v])
        end
        println()
      end
    end
  end

  # additional integrals
  analyze_integrals(analysis_integrals, io, du, u, t, semi)

  return l2_error, linf_error
end

# Print level information only if AMR is enabled
function print_amr_information(callbacks, mesh, solver, cache)

  # Return early if there is nothing to print
  uses_amr(callbacks) || return nothing

  levels = Vector{Int}(undef, nelements(solver, cache))
  min_level = typemax(Int)
  max_level = typemin(Int)
  for element in eachelement(solver, cache)
    current_level = mesh.tree.levels[cache.elements.cell_ids[element]]
    levels[element] = current_level
    min_level = min(min_level, current_level)
    max_level = max(max_level, current_level)
  end

  for level = max_level:-1:min_level+1
    mpi_println(" ├── level $level:    " * @sprintf("% 14d", count(==(level), levels)))
  end
  mpi_println(" └── level $min_level:    " * @sprintf("% 14d", count(==(min_level), levels)))

  return nothing
end

# Print level information only if AMR is enabled
function print_amr_information(callbacks, mesh::P4estMesh, solver, cache)

  # Return early if there is nothing to print
  uses_amr(callbacks) || return nothing

  elements_per_level = zeros(P4EST_MAXLEVEL + 1)

  for tree in unsafe_wrap_sc(p4est_tree_t, mesh.p4est.trees)
    elements_per_level .+= tree.quadrants_per_level
  end

  min_level = findfirst(i -> i > 0, elements_per_level) - 1
  max_level = findlast(i -> i > 0, elements_per_level) - 1

  for level = max_level:-1:min_level+1
    mpi_println(" ├── level $level:    " * @sprintf("% 14d", elements_per_level[level + 1]))
  end
  mpi_println(" └── level $min_level:    " * @sprintf("% 14d", elements_per_level[min_level + 1]))

  return nothing
end

# Iterate over tuples of analysis integrals in a type-stable way using "lispy tuple programming".
function analyze_integrals(analysis_integrals::NTuple{N,Any}, io, du, u, t, semi) where {N}

  # Extract the first analysis integral and process it; keep the remaining to be processed later
  quantity = first(analysis_integrals)
  remaining_quantities = Base.tail(analysis_integrals)

  res = analyze(quantity, du, u, t, semi)
  if mpi_isroot()
    @printf(" %-12s:", pretty_form_utf(quantity))
    @printf("  % 10.8e", res)
    @printf(io, "  % 10.8e", res)
  end
  mpi_println()

  # Recursively call this method with the unprocessed integrals
  analyze_integrals(remaining_quantities, io, du, u, t, semi)
  return nothing
end

# terminate the type-stable iteration over tuples
function analyze_integrals(analysis_integrals::Tuple{}, io, du, u, t, semi)
  nothing
end

# used for error checks and EOC analysis
function (cb::DiscreteCallback{Condition,Affect!})(sol) where {Condition, Affect!<:AnalysisCallback}
  analysis_callback = cb.affect!
  semi = sol.prob.p
  @unpack analyzer = analysis_callback
  cache_analysis = analysis_callback.cache

  l2_error, linf_error = calc_error_norms(sol.u[end], sol.t[end], analyzer, semi, cache_analysis)
  (; l2=l2_error, linf=linf_error)
end

# some common analysis_integrals
# to support another analysis integral, you can overload
# Trixi.analyze, Trixi.pretty_form_utf, Trixi.pretty_form_ascii
function analyze(quantity, du, u, t, semi::AbstractSemidiscretization)
  mesh, equations, solver, cache = mesh_equations_solver_cache(semi)
  analyze(quantity, du, u, t, mesh, equations, solver, cache)
end
function analyze(quantity, du, u, t, mesh, equations, solver, cache)
  integrate(quantity, u, mesh, equations, solver, cache, normalize=true)
end
pretty_form_utf(quantity) = get_name(quantity)
pretty_form_ascii(quantity) = get_name(quantity)

function entropy_timederivative end
pretty_form_utf(::typeof(entropy_timederivative)) = "∑∂S/∂U ⋅ Uₜ"
pretty_form_ascii(::typeof(entropy_timederivative)) = "dsdu_ut"

pretty_form_utf(::typeof(entropy)) = "∑S"

pretty_form_utf(::typeof(energy_total)) = "∑e_total"
pretty_form_ascii(::typeof(energy_total)) = "e_total"

pretty_form_utf(::typeof(energy_kinetic)) = "∑e_kinetic"
pretty_form_ascii(::typeof(energy_kinetic)) = "e_kinetic"

pretty_form_utf(::typeof(energy_kinetic_nondimensional)) = "∑e_kinetic*"
pretty_form_ascii(::typeof(energy_kinetic_nondimensional)) = "e_kinetic*"

pretty_form_utf(::typeof(energy_internal)) = "∑e_internal"
pretty_form_ascii(::typeof(energy_internal)) = "e_internal"

pretty_form_utf(::typeof(energy_magnetic)) = "∑e_magnetic"
pretty_form_ascii(::typeof(energy_magnetic)) = "e_magnetic"

pretty_form_utf(::typeof(cross_helicity)) = "∑v⋅B"
pretty_form_ascii(::typeof(cross_helicity)) = "v_dot_B"

pretty_form_utf(::Val{:l2_divb}) = "L2 ∇⋅B"
pretty_form_ascii(::Val{:l2_divb}) = "l2_divb"

pretty_form_utf(::Val{:linf_divb}) = "L∞ ∇⋅B"
pretty_form_ascii(::Val{:linf_divb}) = "linf_divb"

pretty_form_utf(::typeof(lake_at_rest_error)) = "∑|H₀-(h+b)|"
pretty_form_ascii(::typeof(lake_at_rest_error)) = "|H0-(h+b)|"

# specialized implementations specific to some solvers
include("analysis_dg1d.jl")
include("analysis_dg2d.jl")
include("analysis_dg2d_parallel.jl")
include("analysis_dg3d.jl")
include("analysis_dg3d_parallel.jl")

end # @muladd


text:
8-element Array{Complex{Float64},1}:
      4.0+0.0im
  1.0-2.41421im
      0.0+0.0im
 1.0-0.414214im
      0.0+0.0im
 1.0+0.414214im
      0.0+0.0im
  1.0+2.41421im

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# some commentary

using ..ChemistryFeaturization.AbstractType: AbstractFeaturization
using ..ChemistryFeaturization.FeatureDescriptor:
    AbstractAtomFeatureDescriptor, AbstractPairFeatureDescriptor

struct WeaveFeaturization <: AbstractFeaturization
    atom_features::Vector{<:AbstractAtomFeatureDescriptor}
    pair_features::Vector{<:AbstractPairFeatureDescriptor}
end

function encodable_elements(fzn::WeaveFeaturization)
    # TODO: implement me!
end


text:
using Flux

const PolicyOrSelector = Union{AbstractPolicy, AbstractActionSelector}

mutable struct DQN{Tn<:NeuralNetworkQ, Tp<:PolicyOrSelector, Tf, Tl<:Union{Float32, Float64}} <: AbstractModelFreeLearner
    Q::Tn
    π::Tp
    γ::Float64
    batch_size::Int
    loss_fun::Tf
    loss::Tl
    function DQN(Q::TQ, π::Tp; γ=0.99, batch_size=32, loss_fun=Flux.mse) where {TQ, Tp}
        init_loss = is_using_gpu ? Float32(0.) : Float64(0.)
        new{TQ, Tp, typeof(loss_fun), typeof(init_loss)}(Q, π, γ, batch_size, loss_fun, init_loss)
    end
end

(learner::DQN{<:NeuralNetworkQ, <:AbstractActionSelector})(s) = learner.Q(gpu(s)) |> learner.π
(learner::DQN{<:NeuralNetworkQ, <:AbstractPolicy})(s) = learner.π(s)
(learner::DQN{<:NeuralNetworkQ, <:AbstractPolicy})(s, ::Val{:dist}) = learner.π(s, Val(:dist))

function update!(learner::DQN{<:NeuralNetworkQ, <:AbstractActionSelector}, buffer::CircularSARDBuffer)
    Q, π, γ, batch_size, loss_fun = learner.Q, learner.π, learner.γ, learner.batch_size, learner.loss_fun

    if length(buffer) > batch_size
        (s, a, r, d, s′), _ = sample(buffer, batch_size)
        s, r, d, s′ = gpu(s), gpu(r), gpu(d), gpu(s′)

        q, q′ = Q(s, a), Q(s′, Val(:max))
        G = @. r + γ * q′ * (1 - d)
        loss = loss_fun(G, q)
        learner.loss = loss.data
        update!(Q, loss)
        # update!(π, s, Q(s, Val(:argmax)))  # π isa AbstractPolicy
    end
end

##############################
# Double DQN
##############################

mutable struct DoubleDQN{Ta<:NeuralNetworkQ, Tt<:NeuralNetworkQ, Tp<:AbstractActionSelector, Tf, Tl} <: AbstractModelFreeLearner
    Qₐ::Ta
    Qₜ::Tt
    action_selector::Tp
    update_period::Int
    target_update_period::Int
    training_steps::Int
    batch_size::Int
    loss_fun::Tf
    loss::Tl
    γ::Float64
    function DoubleDQN(Qₐ::Ta, Qₜ::Tt, action_selector::Tp; update_period=4, target_update_period=8000, batch_size=32, loss_fun=Flux.mse, γ=0.99) where {Ta, Tt, Tp}
        init_loss = is_using_gpu ? Float32(0.) : Float64(0.)
        new{Ta, Tt, Tp, typeof(loss_fun), typeof(init_loss)}(Qₐ, Qₜ, action_selector, update_period, target_update_period, 0, batch_size, loss_fun, init_loss, γ)
    end
end

(learner::DoubleDQN)(s) = learner.Qₐ(gpu(s)) |> learner.action_selector

function update!(learner::DoubleDQN, buffer::CircularSARDBuffer)
    Qₐ, Qₜ, γ, batch_size = learner.Qₐ, learner.Qₜ, learner.γ, learner.batch_size
    learner.training_steps += 1

    if length(buffer) > batch_size
        if learner.training_steps % learner.update_period == 0
            (s, a, r, d, s′), _ = sample(buffer, batch_size)
            s, r, d, s′ = gpu(s), gpu(r), gpu(d), gpu(s′)

            q, q′ = Qₐ(s, a), Qₜ(s′, Qₐ(s′, Val(:argmax)))  # note the difference compared to DQN
            G = @. r + γ * q′ * (1 - d)
            loss = learner.loss_fun(G, q)
            learner.loss = loss.data
            update!(Qₐ, loss)
        end

        if learner.training_steps % learner.target_update_period == 0
            Flux.loadparams!(Qₜ.model, Qₐ.ps)
        end
    end
end


text:
#
# national insurance changes from 6th July 2022
# see: https://www.gov.uk/government/publications/spring-statement-2022-factsheet-on-personal-tax/spring-statement-2022-personal-tax-factsheet
#
#
# see fn 5 in that doc on why they say 11,908 and we say 12_570
# for class 4
#
sys.ni.primary_class_1_bands = [123.0, 242.0, 967.0, 9999999999999.9] # the '-1' here is because json can't write inf
sys.ni.class_4_bands = [12_570.0 , 50_270.0, 99999999999999.9 ]

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
# Copyright (c) 2015 Prometheus Research, LLC
#
# Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
# http://www.apache.org/licenses/LICENSE-2.0> or the MIT license
# <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your
# option.  All the files in this project may not be copied, modified,
# or distributed except according to those terms.

#
# Rabbit, a combinator-based query language.
#

__precompile__()

module RBT

include("export.jl")
include("importbase.jl")
include("syntax.jl")
include("immdict.jl")
include("data.jl")
include("type.jl")
include("flow.jl")
include("query.jl")
include("cse.jl")
include("bind.jl")
include("combinator.jl")
include("pgsql.jl")
include("xml.jl")

end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
"""
    struct CoreModelCoupled <: MetabolicModel

The linear model with additional coupling constraints in the form
```
    cₗ ≤ C x ≤ cᵤ
```
"""
mutable struct CoreModelCoupled <: MetabolicModel
    lm::CoreModel
    C::SparseMat
    cl::Vector{Float64}
    cu::Vector{Float64}

    function CoreModelCoupled(lm::MetabolicModel, C::MatType, cl::VecType, cu::VecType)
        length(cu) == length(cl) ||
            throw(DimensionMismatch("`cl` and `cu` need to have the same size"))
        size(C) == (length(cu), n_reactions(lm)) ||
            throw(DimensionMismatch("wrong dimensions of `C`"))

        new(convert(CoreModel, lm), sparse(C), collect(cl), collect(cu))
    end
end

"""
    reactions(a::CoreModelCoupled)

Extract reactions from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
reactions(a::CoreModelCoupled) = reactions(a.lm)

"""
    metabolites(a::CoreModelCoupled)

Extract metabolites from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
metabolites(a::CoreModelCoupled) = metabolites(a.lm)

"""
    stoichiometry(a::CoreModelCoupled)

Extract stoichiometry from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
stoichiometry(a::CoreModelCoupled) = stoichiometry(a.lm)

"""
    bounds(a::CoreModelCoupled)

Extract bounds from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
bounds(a::CoreModelCoupled) = bounds(a.lm)

"""
    balance(a::CoreModelCoupled)

Extract balance from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
balance(a::CoreModelCoupled) = balance(a.lm)

"""
    objective(a::CoreModelCoupled)

Extract objective from [`CoreModelCoupled`](@ref) (uses the internal
[`CoreModel`](@ref)).
"""
objective(a::CoreModelCoupled) = objective(a.lm)

"""
    coupling(a::CoreModelCoupled)::SparseMat

Coupling constraint matrix for a `CoreModelCoupled`.
"""
coupling(a::CoreModelCoupled)::SparseMat = a.C

"""
    n_coupling_constraints(a::CoreModelCoupled)::Int

The number of coupling constraints in a `CoreModelCoupled`.
"""
n_coupling_constraints(a::CoreModelCoupled)::Int = size(a.C, 1)

"""
    coupling_bounds(a::CoreModelCoupled)::Tuple{Vector{Float64},Vector{Float64}}

Coupling bounds for a `CoreModelCoupled`.
"""
coupling_bounds(a::CoreModelCoupled)::Tuple{Vector{Float64},Vector{Float64}} = (a.cl, a.cu)

"""
    reaction_stoichiometry(model::CoreModelCoupled, rid::String)::Dict{String, Float64}

Return the stoichiometry of reaction with ID `rid`.
"""
reaction_stoichiometry(m::CoreModelCoupled, rid::String) = reaction_stoichiometry(m.lm, rid)

"""
    reaction_stoichiometry(model::CoreModelCoupled, ridx)::Dict{String, Float64}

Return the stoichiometry of reaction at index `ridx`.
"""
function reaction_stoichiometry(m::CoreModelCoupled, ridx)::Dict{String,Float64}
    reaction_stoichiometry(m.lm, ridx)
end

"""
    reaction_gene_association_vec(model::CoreModelCoupled)::Vector{Maybe{GeneAssociation}}

Retrieve a vector of gene associations in a [`CoreModelCoupled`](@ref), in the
same order as `reactions(model)`.
"""
reaction_gene_association_vec(model::CoreModelCoupled)::Vector{Maybe{GeneAssociation}} =
    reaction_gene_association_vec(model.lm)

"""
    reaction_gene_association(model::CoreModelCoupled, ridx::Int)::Maybe{GeneAssociation}

Retrieve the [`GeneAssociation`](@ref) from [`CoreModelCoupled`](@ref) by reaction
index.
"""
reaction_gene_association(model::CoreModelCoupled, ridx::Int)::Maybe{GeneAssociation} =
    reaction_gene_association(model.lm, ridx)

"""
    reaction_gene_association(model::CoreModelCoupled, rid::String)::Maybe{GeneAssociation}

Retrieve the [`GeneAssociation`](@ref) from [`CoreModelCoupled`](@ref) by reaction ID.
"""
reaction_gene_association(model::CoreModelCoupled, rid::String)::Maybe{GeneAssociation} =
    reaction_gene_association(model.lm, rid)

"""
    Base.convert(::Type{CoreModelCoupled}, mm::MetabolicModel)

Make a `CoreModelCoupled` out of any compatible model type.
"""
function Base.convert(::Type{CoreModelCoupled}, mm::MetabolicModel)
    if typeof(mm) == CoreModelCoupled
        return mm
    end

    (cl, cu) = coupling_bounds(mm)
    CoreModelCoupled(convert(CoreModel, mm), coupling(mm), cl, cu)
end


text:
using Trilinos
using Test
using MPI

if Test.get_testset_depth() == 0
  MPI.Init()
end

comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))

const n = 20

rowmap0 = Tpetra.Map(n, 0, comm)
@test Tpetra.getIndexBase(rowmap0) == 0

rowmap1 = Tpetra.Map(n, 1, comm)
@test Tpetra.getIndexBase(rowmap1) == 1

if Test.get_testset_depth() == 0
  MPI.Finalize()
end


text:
module AnyMOD

    # XXX enforce use of Julia's own python distribution to avoid interference with local python installations
    using Pkg
    # save current value of environment variable
    if "PYTHON" in keys(ENV)
        envPy = ENV["PYTHON"]
    else
        envPy = ""
    end
    # build python package with Julia distribution
    ENV["PYTHON"]=""
    Pkg.build("PyCall")
    # re-sets environment to former status
    if envPy == ""
        delete!(ENV,"PYTHON")
    else
        ENV["PYTHON"] = envPy
    end

    using Base.Threads, CSV, Dates, LinearAlgebra, Requires
    using MathOptInterface, Reexport, Statistics, PyCall, SparseArrays
    @reexport using DataFrames, JuMP

    pyimport_conda("networkx","networkx")
    pyimport_conda("matplotlib.pyplot","matplotlib")
    pyimport_conda("plotly","plotly")

    include("objects.jl")
    include("tools.jl")
    include("modelCreation.jl")

    include("optModel/exchange.jl")
    include("optModel/objective.jl")
    include("optModel/other.jl")
    include("optModel/tech.jl")

    include("dataHandling/mapping.jl")
    include("dataHandling/parameter.jl")
    include("dataHandling/readIn.jl")
    include("dataHandling/tree.jl")
    include("dataHandling/util.jl")

    export anyModel, initializeModel, createOptModel!, setObjective!
    export reportResults, reportTimeSeries, printObject, printDuals
    export plotTree, plotEnergyFlow, moveNode!

    # XXX define function to print subset of infeasible constraints, if gurobi can be used (has to be installed separately)
    function __init__()
        @require Gurobi="2e9cd046-0924-5485-92f1-d5272153d98b" include("dataHandling/gurobiTools.jl")
    end
end


text:
# # Multiple bodies

#md # ```@meta
#md # CurrentModule = ImmersedLayers
#md # ```

#=
Under the hood, the cache uses the concept of a `Body` (from the `RigidBodyTools.jl` package)
to perform certain
calculations, like normal vectors and surface panel areas, which may
specialize depending on the type of body shape. Note that most immersed layer
operations do not depend on whether there is one or more bodies; rather,
they only depend on the discrete points, and their associated normals and areas.
However, some post-processing operations, like surface integrals, do
depend on distinguishing one body from another. For this reason, the
cache stores points in a `BodyList`, and several operations can exploit this.
=#

using ImmersedLayers
using Plots

#=
For the demonstration, we use the same grid.
=#
Δx = 0.01
Lx = 4.0
xlim = (-Lx/2,Lx/2)
ylim = (-Lx/2,Lx/2)
g = PhysicalGrid(xlim,ylim,Δx)

#=
We will create a 2 x 2 array of circles, each of radius 0.5, centered at $(1,1)$, $(1,-1)$,
$(-1,1)$, $(-1,-1)$.
=#
RadC = 0.5
Δs = 1.4*cellsize(g)
body = Circle(RadC,Δs)

#=
We set up the body list by pushing copies of the same body onto the list.
(We use `deepcopy` to ensure that these are copies, rather than pointers
to the same body.)
=#
bl = BodyList()
push!(bl,deepcopy(body))
push!(bl,deepcopy(body))
push!(bl,deepcopy(body))
push!(bl,deepcopy(body))

#=
Now we move them into position. We also use a `RigidTransform` for each,
which we also assemble into a list. (The `!` is for convenience, using Julia
convention, to remind us that each transform operates in-place on the body.)
=#
t1! = RigidTransform((1.0,1.0),0.0)
t2! = RigidTransform((1.0,-1.0),0.0)
t3! = RigidTransform((-1.0,1.0),0.0)
t4! = RigidTransform((-1.0,-1.0),0.0)
tl! = RigidTransformList([t1!,t2!,t3!,t4!])
nothing #hide

#=
Finally, we apply the transform. We can apply the transform list directly
to the body list:
=#
tl!(bl)

#=
Now we can create the cache, and inspect it by plotting
=#

cache = SurfaceScalarCache(bl,g,scaling=GridScaling)
plot(cache,xlims=(-2,2),ylims=(-2,2))

#=
## Body-by-body calculations
We can now perform operations on data that exploit the division into
distinct bodies. For example, let's compute the integral of $\mathbf{x}\cdot\mathbf{n}$,
for body 3. For any of the bodies, this integral should be approximately equal to the
area enclosed by the body (or volume in 3-d), multiplied by 2 (or 3 in 3-d).
For a circle of radius $1/2$, this area is $\pi/4$, so we expect the result
to be nearly $\pi/2$. We use the `pointwise_dot` operation in
`CartesianGrids.jl` to perform the dot product at each point.
=#
pts = points(cache)
nrm = normals(cache)
V3 = integrate(pointwise_dot(pts,nrm),cache,3)

#=
We can also integrate `VectorData` over individual bodies, and the result
is simply a vector with the integral in each coordinate direction. Let's
demonstrate on another geometric integral, this time of $(\mathbf{x}\cdot\mathbf{x})\mathbf{n}$,
This integral, when divided by the enclosed area (volume) of the body,
is equal to the centroid of the body. Let's demonstrate on body 3, which
we expect to be centered at (-1,1)
=#
Xc = integrate(pointwise_dot(pts,pts)∘nrm,cache,3)/V3

#=
Other operations we can perform body-by-body are [`dot`](@ref) and [`norm`](@ref)
=#

#=
## Copying data body by body
It is common that we will want to assign values to surface data, one body
at a time. For this, we can make use of an extension of the `copyto!` function.
Let's see some examples. Suppose we wish to set the value of a surface
scalar `u` to the $x$ component of the normal vectors for points on body 3, but
leave the values zero for all other bodies. Then we just do the following:
=#
u = zeros_surface(cache)
copyto!(u,nrm.u,cache,3)
nothing #hide

#=
Let's plot the data to verify this worked
=#
plot(u)

#=
It is also possible to use [`copyto!`](@ref) to copy a vector of just the right
size of the subarray associated with the body.
=#


text:
const WP_EST_N = 3
const WP_EST_θ1 = 10 / 180 * pi
const WP_EST_θ2 = 80 / 180 * pi
const F_CUTOFF = 1e-3
const CC_ITER_MAX = 10

function gaplengths(polim::PolarImage, thresh, θ1::Real, θ2::Real)
    checkθ1θ2(θ1,θ2)
    θ1ind = searchsortedfirst(polim.cl.ρ²sort, polim.cl.fθρ(θ1)^2)
    θ2ind = searchsortedlast(polim.cl.ρ²sort, polim.cl.fθρ(θ2)^2)

    threshT = convert(eltype(PolarImage), thresh)
    λ = 0 #gaplength
    gaps = Int[]
    for pixel in polim.imgspiral[θ1ind:θ2ind]
        if pixel < threshT
            λ != 0 && push!(gaps, λ)
            λ = 0
        else
            λ += 1
        end
    end
    return(gaps)
end

# methods for P approach
# Probe probability from formulae 17-18, Chen & Cihlar 1995, IEEE
function probeprob(gaps::Vector{Int})
    length(gaps) == 0 && return([1])
    sortedgaps = sort(gaps)
    len_max = sortedgaps[end]
    P = zeros(Int, len_max)
    for len = 1:len_max
        firstind = searchsortedfirst(sortedgaps, len)
        s = 0
        for λ in sortedgaps[firstind:end]
            s += λ - len
        end
        P[len] = s
    end
    return(P/sum(P)) #normalize probability
end

#Methods for F approach
function reducegaps(gaps::Vector{Int}, Lp, Wp; cutoff=F_CUTOFF)
    # F = theoretical accumulated gap fraction for random foliage (no clumping)
    F = Float64[(1 + Lp*λ/Wp) * exp(-Lp*(1+λ/Wp)) for λ = 0:maximum(gaps)]
    trunc_gaplength = searchsortedfirst(F, cutoff, rev=true)
    sortedgaps = sort(gaps)
    gapsr = sortedgaps[1:searchsortedlast(sortedgaps, trunc_gaplength)]
end
function Fmeasured(gaps::Vector{Int}, pixs)
    gaprange, gapcount = hist(gaps, 0:maximum(gaps))
    blacks = length(pixs) - sum(gaps)
    Fm = 1 - cumsum([blacks, gapcount .* gaprange[2:end]])/length(pixs)
end

# We estimate the projected average leaf width once per image, because
# it's not stable for small zenith angles θ and "the calculation of
# [clumping index] Ω with a gap-size distribution is not very sensitive to the
# choice of the elementh width within a reasonable range" (Chen,Chilar 1995,
# applied optics, section 4.B).
# Wp estimation according to Leblanc et al, 2005, formula (7).
function Wp_estimation(polim::PolarImage, thresh;
    θ1=WP_EST_θ1, θ2=WP_EST_θ2, Wp_est_n=WP_EST_N)

    checkθ1θ2(θ1,θ2)
    gaps = gaplengths(polim, thresh, θ1, θ2)
    Wp =  Wp_estimation(probeprob(gaps); n=Wp_est_n)
end

function Wp_estimation(P; n=WP_EST_N)
    regrcoef = linreg([1.:n], log(P[1:n]))
    Wp = - regrcoef[1] / abs(regrcoef[2])
end

function chencihlar(polim::PolarImage, thresh, θ1::Real, θ2::Real; kwargs...)
    checkθ1θ2(θ1,θ2)
    # Wp estimation not stable for small zenith angles, so estimate Wp from full image
    if θ1 > WP_EST_θ1
        #TODO add Wp_est_n keyword argument
        Wp = Wp_estimation(probeprob(gaplengths(polim, thresh, θ1, θ2)))
    else
        Wp = Wp_estimation(polim, thresh)
    end
    chencihlarF(Wp, polim, thresh, θ1, θ2; kwargs...)
end

# P method Chen Cihlar, requires Wp input
function chencihlarF(Wp::Float64, polim::PolarImage, thresh, θ1::Real, θ2::Real;
    cutoff=F_CUTOFF, Wp_est_n=WP_EST_N, itermax=CC_ITER_MAX)

    checkθ1θ2(θ1,θ2)
    gaps = gaplengths(polim, thresh, θ1, θ2)
    length(gaps) == 0 && return(Ω = 1)
    pixs = pixels(polim, θ1, θ2)

    # Lp initiation
    Fm = Fmeasured(gaps, pixs)
    Lp = -log(Fm[1])

    # loop Lp
    Lp_prev = Lp
    counter = 0
    Fmr = Float64[] # initiation required before loop
    while counter < itermax
        counter += 1
        gaps = reducegaps(gaps, Lp, Wp)
        Fmr = Fmeasured(gaps, pixs)
        Lp = -log(Fmr[1])
        abs(1-Lp/Lp_prev) < 0.01 && break
        Lp_prev = Lp
    end
    counter == itermax && @warn("chencihlar max iteration reached")
    Ω = log(Fm[1]) / log(Fmr[1]) * (1-Fmr[1])/(1-Fm[1])
end

# P-method from Chen & Cihlar, IEEE trans. geo rem. sens. 1995
# Gives low clumping factors.
function chencihlarP(polim::PolarImage, thresh, θ1::Real, θ2::Real;n=WP_EST_N)
    checkθ1θ2(θ1,θ2)
    gaps = gaplengths(polim, thresh, θ1, θ2)
    Pl = probeprob(gaps)
    length(Pl) < 2 && return(1)
    maxP = length(Pl)
    n = min(n, maxP)
    Leθ = -linreg([1.:n], log(Pl[1:n]))[1]

    Pstart = ifloor(maxP/4)
    Pstop = ifloor(maxP/3)

    Lcθ = -linreg(float([Pstart:Pstop]), log(Pl[Pstart:Pstop]))[1]
    α = 0
    LEp = log((1+α)*Lcθ*exp(-Lcθ)/(sqrt(2(1+α)exp(-(Leθ+Lcθ))-(1+2α)exp(-2Lcθ))-exp(-Lcθ)))
    α = Lcθ*exp(-LEp)/3 #higher order correction
    LEp = log((1+α)*Lcθ*exp(-Lcθ)/(sqrt(2(1+α)exp(-(Leθ+Lcθ))-(1+2α)exp(-2Lcθ))-exp(-Lcθ)))
    Ω = Leθ / (LEp*Lcθ)
end


text:
# parsemat.jl
#
# parses impedance matrix output file from FastFieldSolvers
export parsefasthenrymat

"""
    ParseFastHenryMatResult

Result of parsing FastHenry .mat file.

**Fields**

- `portnames`    -- array of port names, index is row number in impedance matrix
- `frequencies`  -- frequencies at which impedance matrix is computed
- `impedance`    -- impedance matrix at each frequency. impedance[row, col, frequency]

"""
mutable struct ParseFastHenryMatResult
  "array of port names, index is row number in impedance matrix"
  portnames :: Array{AbstractString,1}
  "list of frequencies at which impedance matrix is computed"
  frequencies :: Array{Float64,1}
  """
  impedance matrix at each frequency.  <br>
  impedance [impedance matrix row : impedance matrix column : frequency]
  """
  impedance :: Array{ComplexF64,3}
end

function parsefasthenrymat(io::IO)
  findportname = r"Row (\d+):.*?(port name: (\w+)|$)"m
  line = readline(io)
  m = match(findportname,line)
  if m==nothing
    throw(Meta.ParseError("first line did not start with \"Row\""))
  end
  numberofrows = parse(Int,m.captures[1])
  portnames = Array{String}(undef,numberofrows)
  portnames[numberofrows] = m.captures[3] == nothing ? "" : m.captures[3]
  for i in numberofrows-1:-1:1
    line = readline(io)
    m = match(findportname,line)
    if m==nothing
      rowstring = @sprintf("%d",i)
      throw(Meta.ParseError("did not find rows line for Row="*rowstring))
    end
    portnames[i] = m.captures[3] == nothing ? "" : m.captures[3]
  end
  frequencies = 0
  buffer = IOBuffer()
  impeadancematrix = Array{ComplexF64,2}(undef,numberofrows,numberofrows)
  while ~eof(io)
    line = readline(io)
    frequencies+=1
    m = match(r"Impedance matrix for frequency = ([\d+-eE.]+) ",line)
    if m==nothing
      throw(Meta.ParseError("did not find \"Impedance matrix for frequency = \""))
    end
    frequency = parse(Float64,m.captures[1])
    write(buffer,frequency)
    # read impedance matrix here
    for row in 1:numberofrows
      line = readline(io)
      position = 1
      for column in 1:numberofrows # impedance matrix is square
        m = match(r"([\d+-.eE]+) *([\d+-.eE]+)j",line,position)
        if m==nothing
          throw(Meta.ParseError("error parsing impedance matrix"))
        end
        r = parse(Float64,m.captures[1])
        i = parse(Float64,m.captures[2])
        impeadancematrix[row,column]=ComplexF64(r,i)
        position = m.offset + length(m.match)
      end
    end
    write(buffer,impeadancematrix)
  end
  resultmatrix = Array{ComplexF64,3}(undef,numberofrows,numberofrows,frequencies)
  frequencymatrix = Array{Float64}(undef,frequencies)
  seekstart(buffer)
  for i in 1:frequencies
    frequency = read(buffer,Float64)
    frequencymatrix[i] = frequency
    read!(buffer,impeadancematrix)
    resultmatrix[:,:,i] = impeadancematrix
  end
  ParseFastHenryMatResult(portnames,frequencymatrix,resultmatrix)
end

function parsefasthenrymat(filename::String)
  io = open(filename,"r")
  result = parsefasthenrymat(io)
  close(io)
  return result
end

"""
    parsefasthenrymat(io::IO)
    parsefasthenrymat(filename::AbstractString)

parses the .mat output file from FastHenry

returns type `ParseFastHenryResult`
"""
parsefasthenrymat


text:
include("../perfutil.jl")

include("level1.jl")
include("level2.jl")
include("level3.jl")

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using BinningAnalysis
using Test, Statistics, Random

@testset "All Tests" begin

    @testset "Logarithmic Binning" begin
        include("logbinning.jl")
    end

    @testset "Full Binning" begin
        include("fullbinning.jl")
    end

    @testset "Error Propagator" begin
        include("errorpropagator.jl")
    end

    @testset "Incremental Binning" begin
        include("incrementbinner.jl")
    end

    @testset "Generic functions" begin
        x = [0.0561823, 0.846613, 0.813439, 0.357134, 0.157445, 0.103298, 0.948842, 0.629425, 0.290206, 0.00695332, 0.869828, 0.949165, 0.897995, 0.916239, 0.457564, 0.349827, 0.398683, 0.264218, 0.72754, 0.934315, 0.666448, 0.134813, 0.364933, 0.829088, 0.256443, 0.595029, 0.172097, 0.241686, 0.489935, 0.239663, 0.391291, 0.00751015, 0.138935, 0.0569876, 0.571786, 0.694996, 0.798602, 0.923308, 0.73978, 0.414774, 0.835145, 0.731303, 0.271647, 0.707796, 0.00348624, 0.0905812, 0.316176, 0.921054, 0.131037, 0.599667, 0.805071, 0.440813, 0.086516, 0.363658, 0.476161, 0.931257, 0.28974, 0.78717, 0.60822, 0.144024, 0.214432, 0.061922, 0.626495, 0.512072, 0.758078, 0.840485, 0.242576, 0.147441, 0.599222, 0.993569, 0.0365044, 0.0983033, 0.713144, 0.422394, 0.480044, 0.968745, 0.518475, 0.431319, 0.4432, 0.526007, 0.612975, 0.468387, 0.262145, 0.888011, 0.105744, 0.325821, 0.769525, 0.289073, 0.336083, 0.443037, 0.489698, 0.141654, 0.915284, 0.319068, 0.341001, 0.704346, 0.0794996, 0.0412352, 0.70016, 0.0195158]
        @test std_error(x) ≈ std_error(x; method=:log)
        @test std_error(x; method=:log) ≈ 0.03384198827323159
        @test std_error(x; method=:full) ≈ 0.030761063167290947
        @test std_error(x; method=:jackknife) ≈ 0.029551526503551463
        @test std_error(x; method=:error_propagator) ≈ 0.03384198827323159
        @test tau(x) ≈ tau(x; method=:log)
        @test tau(x; method=:log) ≈ 0.15572524869061377
        @test tau(x; method=:full) ≈ 0.04176737474771308

        x = x .+ 1im
        @test std_error(x)≈ std_error(x; method=:log)
        @test std_error(x; method=:log) ≈ 0.0338419882732316
        @test std_error(x; method=:full) ≈ 0.03076106316729094
        @test std_error(x; method=:jackknife) ≈ 0.029551526503551445
        @test std_error(x; method=:error_propagator) ≈ 0.03384198827323173
        @test tau(x) ≈ tau(x; method=:log)
        @test tau(x; method=:log) ≈ 0.15572524869061377
        @test tau(x; method=:full) ≈ 0.04176737474771308

        x = Array{Float64,2}[[0.311954 0.476706 0.314101; 0.420978 0.478085 0.0194284], [0.491386 0.493583 0.295477; 0.136896 0.1634 0.84641], [0.318536 0.705903 0.121377; 0.764174 0.240484 0.25894], [0.286923 0.172458 0.392881; 0.124348 0.140628 0.730131], [0.287451 0.221914 0.382938; 0.29568 0.249575 0.87685]]
        @test std_error(x) ≈ std_error(x; method=:log)
        @test std_error(x, method=:log) ≈ [0.03856212859659066 0.09765048982846934 0.04879645809318544; 0.11744109653881814 0.05978167652935139 0.1722342286233907]
        @test std_error(x, method=:full) ≈ [0.03856212859659072 0.09765048982846936 0.048796458093185405; 0.11744109653881814 0.05978167652935137 0.17223422862339066]
        @test std_error(x, method=:error_propagator) ≈ [0.03856212859659066 0.09765048982846934 0.04879645809318544; 0.11744109653881814 0.05978167652935139 0.1722342286233907]
        @test tau(x) ≈ tau(x; method=:log)
        @test tau(x; method=:log) ≈ [0.0 0.0 0.0; 0.0 0.0 0.0]
        @test isapprox(tau(x; method=:full), [0.0 -5.55112e-17 -5.55112e-17; -1.11022e-16 0.0 0.0], atol=1e-16)

        x = ["this", "should", "error"]
        @test_throws ErrorException std_error(x, method=:jackknife)
        @test_throws ArgumentError std_error(x, method=:whatever)
    end

    @testset "Jackknife" begin
        Random.seed!(123)
        g(x1, x2) = x1^2 - x2

        # Real
        ts = rand(Float64, 1000)
        ts2 = rand(Float64, 1000)

        @test isapprox.(
            jackknife_full(identity, ts),
            (0.5032466005013134, 1.1091128016005314e-13, 0.008897060002353763),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(identity, 1 ./ ts),
            (7.97707906607684, -5.323741447682551e-12, 2.522127468106169),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(g, ts, ts2.^2),
            (-0.05969375347674344, 7.915767673482427e-5, 0.013117185693279606),
            atol = 1e-12
        ) |> all

        # check consistency with Julia's Statistics
        @test jackknife(identity, ts)[2] ≈ std(ts)/sqrt(length(ts))

        # Complex
        ts = rand(ComplexF64, 1000)
        ts2 = rand(ComplexF64, 1000)

        @test isapprox.(
            jackknife_full(identity, ts),
            (0.5026931517594448 + 0.49854565281725627im, 4.4364512064021255e-13 - 3.327338404801594e-13im, 0.012946143510730868),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(identity, 1 ./ ts),
            (1.1496877962692906 - 1.0838709685240246im, 1.3309353619206377e-12 - 6.654676809603188e-13im, 0.05676635162628594),
            atol = 1e-12
        ) |> all
        @test isapprox.(
            jackknife_full(g, ts, ts2.^2),
            (0.02004836872405491 - 0.01825289014938747im, -3.1784715169222433e-6 + 5.949834723031899e-7im, 0.026446578773003968),
            atol = 1e-12
        ) |> all

        # check consistency with Julia's Statistics
        @test jackknife(identity, ts)[2] ≈ std(ts)/sqrt(length(ts))
    end
end


text:
using RandomMatrices
using CUDA

function tensor_init_pauli(; enable_cuda = false)
    tn = Float64[1, 0, 0, 1] / sqrt(2)
    enable_cuda && return CuArray(tn)
    return tn
end
function tensor_measure_pauli(; enable_cuda = false)
    tn = (1/sqrt(2))*Float64[1 1; 0 0; 0 0; 1 -1]
    enable_cuda && return CuArray(tn)
    return tn
end
function tensor_single_pauli(g_prev = :null; haar = false, enable_cuda = false)
    I = [1 0; 0 1]
    X = [0 1; 1 0]
    Y = [0 -im; im 0]
    Z = [1 0; 0 -1]
    paulis = (I, X, Y, Z)
    if haar
        M = rand(Haar(2), 2)
        g_new = :null
    else    # Haar random matrix
        g_prev === :null && (g_new = rand((:sx, :sy, :sw)))
        g_prev === :sx && (g_new = rand((:sy, :sw)))
        g_prev === :sy && (g_new = rand((:sx, :sw)))
        g_prev === :sw && (g_new = rand((:sx, :sy)))

        g_new === :sx && (M = sqrt(X))
        g_new === :sy && (M = sqrt(Y))
        g_new === :sw && (M = sqrt(sqrt(1/2)*(X+Y)))
        M = [1 0; 0 exp(rand()*2π*im)]*M*[1 0; 0 exp(rand()*2π*im)]
    end
    ts = zeros(Float64, 4, 4)
    for i = 1:4, j = 1:4
        ts[i, j] = real(tr(M*paulis[i]*M'*paulis[j]/2))
    end
    enable_cuda && return CuArray(ts), g_new
    return ts, g_new
end
function tensor_noise_pauli(; enable_cuda = false)
    tn = zeros(Float64, 4)
    tn[1] = 1
    enable_cuda && return CuArray(tn)
    return tn
end
function tensor_fsim_pauli(; enable_cuda = false)
    fsim_gate = zeros(ComplexF64, 4, 4)
    fsim_gate[1, 1] = 1
    fsim_gate[2, 3] = -im
    fsim_gate[3, 2] = -im
    fsim_gate[4, 4] = exp(-pi/6*im)
    tf = zeros(Float64, 4, 4, 4, 4)
    I = [1 0; 0 1]
    X = [0 1; 1 0]
    Y = [0 -im; im 0]
    Z = [1 0; 0 -1]
    paulis = (I, X, Y, Z)
    for a = 1:4, b = 1:4, c = 1:4, d = 1:4
        tf[a, b, c, d] = real(tr(fsim_gate * kron(paulis[a], paulis[b]) * fsim_gate' * kron(paulis[c], paulis[d])) / 4)
    end
    enable_cuda && return CuArray(tf)
    return tf
end

"""
    to_ein_code_pauli(layout, [s; haar = false, enable_cuda = false])

Generate the tensor network description for `layout` on a subset `s`. If the is a FSim gate out side `s`,
it will be replaced by an identity gate.
"""
function to_ein_code_pauli(layout::RQCLayout{VT, PT, SCGate}, s = collect(vertices(layout)); haar = false, enable_cuda = false) where {VT, PT}
    layout = deepcopy(layout)
    update_id!(layout)
    vs = sort!(intersect(s, collect(vertices(layout))))
    nbits = length(vs)
    ids_size = Dict{Int, Int}(i => 4 for i in 1:nbits)
    ids_loc = Dict(i => (vs[i], 0) for i in 1:nbits)
    open_ids_dict = Dict(vs[i] => i for i = 1:nbits)
    tensor_ids = [[id] for id in 1:nbits]
    max_id = nbits
    tensors = Any[]
    for v in vs
        if gates(layout, v)[1] === Noise
            push!(tensors, sqrt(1/2)*tensor_noise_pauli(enable_cuda = enable_cuda))
            gates(layout, v)[1] = Id
        else
            push!(tensors, tensor_init_pauli(enable_cuda = enable_cuda))
        end
    end
    last_single = Dict(vs[i] => :null for i = 1:nbits)
    D = length(gates(layout, 1))
    for d = 1:D
        for v in vs
            if gates(layout, v)[d] === Single
                ts, last_single[v] = tensor_single_pauli(last_single[v]; haar = haar, enable_cuda = enable_cuda)
                max_id += 1
                new_id_v = max_id
                push!(tensors, ts)
                push!(tensor_ids, [open_ids_dict[v], new_id_v])
                ids_size[new_id_v] = 4
                ids_loc[new_id_v] = (v, d)
                open_ids_dict[v] = new_id_v
            elseif gates(layout, v)[d] === Noise && d < D
                max_id += 1
                new_id_v = max_id
                push!(tensors, tensor_noise_pauli(enable_cuda = enable_cuda), tensor_noise_pauli(enable_cuda = enable_cuda))
                push!(tensor_ids, [open_ids_dict[v]], [new_id_v])
                ids_size[new_id_v] = 4
                ids_loc[new_id_v] = (v, d)
                open_ids_dict[v] = new_id_v
            elseif gates(layout, v)[d] === FSim
                c = depth_to_cycle(d)
                u = partner(layout, v, c)
                if u < v && u in vs
                    max_id += 2
                    new_id_u = max_id - 1
                    new_id_v = max_id
                    push!(tensors, tensor_fsim_pauli(enable_cuda = enable_cuda))
                    push!(tensor_ids, [open_ids_dict[u], open_ids_dict[v], new_id_u, new_id_v])
                    ids_size[new_id_u] = 4
                    ids_size[new_id_v] = 4
                    ids_loc[new_id_u] = (u, d)
                    ids_loc[new_id_v] = (v, d)
                    open_ids_dict[u] = new_id_u
                    open_ids_dict[v] = new_id_v
                end
            end
        end
    end
    for v in vs
        if gates(layout, v)[end] === Noise
            push!(tensors, sqrt(2)*tensor_noise_pauli(enable_cuda = enable_cuda))
            push!(tensor_ids, [open_ids_dict[v]])
            open_ids_dict[v] = 0
        else
            max_id += 1
            new_id_v = max_id
            push!(tensors, tensor_measure_pauli(enable_cuda = enable_cuda))
            push!(tensor_ids, [open_ids_dict[v], new_id_v])
            ids_size[new_id_v] = 2
            ids_loc[new_id_v] = (v, D)
            open_ids_dict[v] = new_id_v
        end
    end

    open_ids = Int[]
    for v in vs
        open_ids_dict[v] > 0 && push!(open_ids, open_ids_dict[v])
    end
    # return tensor_ids, open_ids
    return EinCode(tensor_ids, open_ids), tensors, ids_size, ids_loc
end