text:
##########
##########
# Definition of functions used to look at species rewiring in relevant paper
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2021-04-29
##########
##########

# set up =======================================================================

# set packages
using DataFrames, Distributions, Random, StatsBase, Documenter, Colors
using LightGraphs, GraphPlot, GraphRecipes, Plots

# trait-related functions ======================================================

"""
    assign_trait_values(organism_type::String, organism_id::Int)

Takes in the type of organism and it's unique id and uses the percentage of its
own body size as defined, to return a list of the various traits for that
organism to be passed to the next list.

"""
function assign_trait_values(
    organism_type::String,
    organism_id::Int,
    lower_prey_limit::Float64 = 0.001,
    upper_prey_limit::Float64 = 0.65
    )
    # determine producer or consumer
    if organism_type == "producer"
        # scale for producers should be 0.34 (Brose 2019 paper on body size)
        body_size = rand(Exponential(0.34), 1)
        # set prey limits - empty for producers
        prey_limits = Array([NaN, NaN])
    else # everything here down is for consumers
        # scale for consumers should be 800 (Brose 2019 paper on body size)
        body_size = rand(Exponential(20), 1)
        # set prey limits - empty for producers
        lower_limit = body_size*lower_prey_limit
        upper_limit = body_size*upper_prey_limit
        prey_limits = [lower_limit[1], upper_limit[1]]

    end
    # habitat values get set by same process for both consumers and producers

    # randomly select habitat midpoint
    habitat_midpoint = round(rand(Uniform(0,1), 1)[1], digits = 3)
    # define vector of the midpoint and upper/lower habitat bounds
    habitat = [habitat_midpoint, 0, 0]
    if habitat_midpoint >= 0.4 # for lower bound
        habitat[2] = round((habitat_midpoint - 0.4), digits = 3)
    else
        habitat[2] = 0
    end
    if habitat_midpoint <= 0.6 # for upper bound
        habitat[3] = round((habitat_midpoint + 0.4), digits = 2)
    else
        habitat[3] = 1
    end
    #nutritional value
    nutritional_value = rand(Beta(5, 2))

    #habitat = [habitat_midpoint, 0, 1]
    trait_list = [organism_id, organism_type, round(body_size[1], digits = 2),
                    round(prey_limits[1], digits = 2),
                    round(prey_limits[2], digits = 2), habitat[1],
                    round(habitat[2], digits = 2),
                    round(habitat[3], digits = 2),
                    round(nutritional_value, digits = 2)]

    return trait_list
end

"""
    define_species_list(n::Int, shape::String)

Takes in the number of species and a shape for the web and returns an array
with each element being the trait list for each species in the web

# examples & checks
```jldoctest
# make sure ArgumentError works
julia> define_species_list(50, "topheavy")
ERROR: ArgumentError: shape is not recognized: can only be 'top-heavy',
'bottom-heavy', or 'uniform'

```
"""
function define_species_list(n::Int, richness_shape::String)
    # make sure shape is of the right kind
    if richness_shape ∉ ["top-heavy", "bottom-heavy", "uniform"]
        ArgumentError("richness_shape is not recognized: can only be 'top-heavy',
        'bottom-heavy', or 'uniform'")
    end
    # define number of species in the array with predefined types of the sublist
    species_list = Array{Any}(undef, 50, 9)
    # set items in organism type list to be selected from
    types = ["producer", "consumer"]
    # if statements to determine the weights of these
    if richness_shape == "top-heavy"
        weights = [1, 5]
    elseif richness_shape == "bottom-heavy"
        weights = [5,1]
    else richness_shape == "uniform"
        weights = [1,1]
    end
    # use while loop to populate list
    i = 1
    while i <= n
        organism_type = sample(types, Weights(weights))
        species_list[i, :] = assign_trait_values(organism_type, i)

        i += 1 #increment
    end
    return species_list
end

"""
    get_species(species_list, i)

Finds the species list entry according to species id, NOT simply location in the
array
"""
function get_species(species_list::Matrix, id::Int)

    # find focal species in array
    matrix_sp = species_list[species_list[:, organism_id_pos] .== id, :]
    focal_sp = []
    focal_sp = append!(focal_sp, matrix_sp)

    return focal_sp
end

"""
    determine_interaction(predator_id::Int, prey_id::Int, species_list::Array)

Assigns whether or not a given predator & prey have an interaction according to
the stated rules of the algorithm

"""
function determine_interaction(
    predator_id::Int,
    prey_id::Int,
    species_list::Array
    )
    # extract the values for habitat, size, and predation limits for ease
    predator = get_species(species_list, predator_id)
    prey = get_species(species_list, prey_id)

    pred_hab = predator[habitat_midpoint_pos:habitat_upper_limit_pos]
    prey_hab = prey[habitat_midpoint_pos:habitat_upper_limit_pos]
    prey_size = prey[body_size_pos]
    pred_min = predator[lower_prey_limit_pos]
    pred_max = predator[upper_prey_limit_pos]

    if (predator_id == prey_id) || (predator[organism_type_pos] == "Producer")
        value = 0 # assuming no canibalism and that producers don't predate
    elseif (((prey_hab[2] >= pred_hab[2] && prey_hab[2] <= pred_hab[3]) ||
        (prey_hab[3] >= pred_hab[2] && prey_hab[3] <= pred_hab[3])) &&
        (prey_size[1] <= pred_max && prey_size[1] >= pred_min))
        value = 1 # assign 1 if within body size and habitat range
    else
        value = 0
    end

    return value
end

"""
    biomass_interaction(predator_id::Int, prey_id::Int, species_list::Array)

Assigns whether or not a given predator & prey have an interaction according the
same stated rules as determine_interaction(), but now including whether either
the predator or prey are extinct (i.e. biomass = 0.0)

"""
function biomass_interaction(
    predator_id::Int,
    prey_id::Int,
    species_list::Array
    )
    # extract the values for habitat, size, and predation limits for ease
    predator = get_species(species_list, predator_id)
    prey = get_species(species_list, prey_id)

    pred_hab = predator[habitat_midpoint_pos:habitat_upper_limit_pos]
    prey_hab = prey[habitat_midpoint_pos:habitat_upper_limit_pos]
    prey_size = prey[body_size_pos]
    pred_min = predator[lower_prey_limit_pos]
    pred_max = predator[upper_prey_limit_pos]

    if (predator_id == prey_id) || (predator[organism_type_pos] == "Producer")
        value = 0 # assuming no canibalism and that producers don't predate
    elseif predator[biomass_pos] == 0.0
        value = 0
    elseif prey[biomass_pos] == 0.0
        value = 0
    elseif (((prey_hab[2] >= pred_hab[2] && prey_hab[2] <= pred_hab[3]) ||
        (prey_hab[3] >= pred_hab[2] && prey_hab[3] <= pred_hab[3])) &&
        (prey_size[1] <= pred_max && prey_size[1] >= pred_min))
        value = 1 # assign 1 if within body size and habitat range
    else
        value = 0
    end

    return value
end

"""
    create_interaction_matrix(species_list::Array)

Takes in the species list and returns an interaction matrix that will form the
basis of the networks to then be created.

"""
function create_interaction_matrix(species_list::Array)
    interaction_matrix = fill(0, (size(species_list, 1), size(species_list, 1)))
    for i in 1:size(interaction_matrix)[1] #rows
        for j in 1:size(interaction_matrix)[2] #cols
            interaction_matrix[i,j] = determine_interaction(j, i, species_list)
        end
    end
    return interaction_matrix
end

"""
    biomass_interaction_matrix(species_list::Array)

Takes in the species list and returns an interaction matrix that will allow for
inspection of networks after biomass simulations have occured.

"""
function biomass_interaction_matrix(species_list::Array)
    interaction_matrix = fill(0, (size(species_list, 1), size(species_list, 1)))
    for i in 1:size(interaction_matrix)[1] #rows
        for j in 1:size(interaction_matrix)[2] #cols
            interaction_matrix[i,j] = biomass_interaction(j, i, species_list)
        end
    end
    return interaction_matrix
end

"""
    check_connectedness(species_list)

Check the connectedness of the interaction_matrix formed by a specific
species list draw.
"""
function check_connectedness(
    species_list::Array,
    )

    interaction_matrix = create_interaction_matrix(species_list)
    # get indices of the consumers to check if they are consuming something
    consumers = species_list[species_list[:,organism_type_pos] .== "consumer",
                                organism_id_pos]
    producers = species_list[species_list[:,organism_type_pos] .== "producer",
                                organism_id_pos]
    producers = species_list[species_list[:,organism_type_pos] .== "g",
                                organism_id_pos]
    # make sure all consumers are eating at least one other organism (i.e. the
    # sum of their column is > 0)
    consumer_matrix = interaction_matrix[:,consumers]
    if size(consumer_matrix, 1) > 0
        consumer_sums = minimum(sum(consumer_matrix, dims = 1))
    else
        return "fail" # this makes sure there are no all-producer networks
    end
    # make sure all producers are being consumeed by least one other organism
    # (i.e. the sum of their column is > 0)
    producer_matrix = interaction_matrix[producers,:]
    if size(producer_matrix, 1) > 0
        producer_sums =  minimum(sum(producer_matrix, dims = 2))
    else
        producer_sums = 1 # just make it non-zero so it will pass the test
    end

    if consumer_sums != 0 && producer_sums != 0
        return "pass"
    else
        return "fail"
    end
end

"""
    network_draw(n::Int64, richness_shape::String)

This function takes in the number of species and the richness shape and draws
an initial network and makes the checks to see if it's valid - if it is
connected, and then resamples the connectance to prepare it to have the network
metrics sampled.
"""
function network_draw(
    n::Int64,
    richness::String
    )

    # draw species list
    species_list = define_species_list(n, richness)
    interaction_matrix = create_interaction_matrix(species_list)
    # check connectedness, and go through again if need be
    if check_connectedness(species_list) == "pass"
        #print("first time success ")
        return species_list, interaction_matrix
    else
        #print("first time fail ")
        # go through while loop to re-draw until connected
        connect_acomplish = false
        connected_counter = 0
        while connect_acomplish == false
            connected_counter += 1
            # draw species list
            species_list = define_species_list(n, richness)
            interaction_matrix = create_interaction_matrix(species_list)
            if check_connectedness(species_list) == "fail"
                continue
            elseif check_connectedness(species_list) == "pass"
                connect_acomplish = true
                #print("success after ", connected_counter, " additional draws.")
                return species_list, interaction_matrix
            end
        end
    end
end

"""
    get_thirds(sequence::Array, num::Int64)

Takes in sequence and an integer and returns the sequence split into `num`
groups as an array of arrays.
"""
function get_thirds(sequence::Array, num::Int64)

    k, m = divrem(length(sequence), num)
    thirds = []
    for i in 0:(num-1)
        append!(thirds, [sequence[(i*k + min(i, m) +1):((i+1)*k + min(i+1, m))]])
    end
    thirds
    return thirds
end

"""
    initialize_biomass(species_list::Array, biomass_shape::String="pyramid")

Initializes the amount of biomass for each species in the species list according
to the shape of the web (i.e. pyramid, uniform, inverted)

"""
function initialize_biomass(
    species_list::Array,
    biomass_shape::String
    )

    # get all body body sizes
    body_sizes = zeros(0)
    for i in 1:size(species_list,1)
        append!(body_sizes, species_list[i,body_size_pos])
    end
    # make three equal groups of the body sizes
    sorted_body_sizes = sort(body_sizes)
    thirds = get_thirds(sorted_body_sizes, 3)
    small = [thirds[1][1],thirds[1][end]]
    medium = [thirds[2][1],thirds[2][end]]
    large = [thirds[3][1],thirds[3][end]]

    if biomass_shape == "pyramid"
        small_biomass = 12.5
        medium_biomass = 5
        large_biomass = 1.25
    elseif biomass_shape == "inverted"
        small_biomass = 1.25
        medium_biomass = 5
        large_biomass = 12.5
    else biomass_shape == "uniform"
        small_biomass = 6.2333333
        medium_biomass = 6.2333333
        large_biomass = 6.2333333
    end

    # get mean biomass for each size bin
    num_small_sp = length(thirds[1])
    num_medium_sp = length(thirds[2])
    num_large_sp = length(thirds[3])

    large_mean = large_biomass/num_large_sp
    medium_mean = medium_biomass/num_medium_sp
    small_mean = small_biomass/num_small_sp

    # get random biomass sample and append it to the end of each species trait info
    temp = zeros(50)
    for i in 1:size(species_list,1)
        if species_list[i,3] ≤ small[2]
            temp[i] = rand(TruncatedNormal(small_mean, 0.09, 0.05, 1))
        elseif species_list[i,3] >= large[1]
            temp[i] = rand(TruncatedNormal(large_mean, 0.09, 0.05, 1))
        else
            temp[i] = rand(TruncatedNormal(medium_mean, 0.09, 0.05, 1))
        end
    end
    species_list = hcat(species_list, temp)

    return species_list
end

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
"""
    spy(A; kw...)

# Description

Plots the sparsity pattern for the given matrix `A`. This means
that a scatterplot that resembles the matrix is drawn, in which
only the pixel for non-zero elements of the matrix are set.

If the parameters `width` and `height` are not explicitly
specified, then the function will attempt to preserve the aspect
ratio of the matrix, while also attempting to fit the resulting
plot withing the bounding box specified by `maxwidth` and `maxheight`.

# Usage

    spy(A; $(keywords((maxwidth = 0, maxheight = 0, title = "Sparsity Pattern", zeros = false); add=(:fix_ar, :canvas))))

# Arguments

$(arguments(
    (
        A = "matrix of interest for which non-zero elements should be drawn",
        maxwidth = "maximum number of characters per row that should be used for plotting",
        maxheight = "maximum number of character rows that should be used for plotting",
        width = "exact number of characters per row that should be used for plotting (`0` stands for automatic)",
        height = "exact number of character rows that should be used for plotting (`0` stands for automatic)",
        show_zeros = "show zeros pattern instead of default nonzeros",
    ); add=(:fix_ar, :canvas),
))

# Returns

A plot object of type `Plot{T<:Canvas}`.

# Author(s)

- Dominique Orban (github.com/dpo)
- Christof Stocker (github.com/Evizero)
- Jake Bolewski (github.com/jakebolewski)

# Examples

```julia-repl
julia> spy(sprandn(50, 120, .05))
                           Sparsity Pattern
      ┌────────────────────────────────────────────────────────────┐
    1 │⢀⠀⠒⠀⢄⠂⠀⡀⠀⠀⠀⠰⠀⠐⠀⠀⠀⠀⠂⠀⠀⠀⠠⠰⠀⠀⠀⡀⠄⠀⠠⠀⠀⠀⠀⠀⡀⠀⠀⠀⡀⠐⠀⠐⠀⠀⡀⠀⠀⠀⠀⠀⠀⠂⡀⠀⡀⠠⠀⠀│ > 0
      │⠀⠀⠀⠀⠀⢀⢀⠀⠀⠀⠀⠐⠄⠀⢀⠂⠀⠠⠀⠉⠀⠀⠀⠀⢐⠀⠀⠀⠀⠀⠠⠐⠄⠄⠀⠀⠀⠐⠐⠈⢁⠀⠀⠀⠀⠀⠀⠀⠠⠀⠐⡐⠀⠄⠀⠀⠈⠀⠀⠠│ < 0
      │⡀⠀⠀⠀⡀⠊⠀⠀⠀⢀⢀⠀⠄⠀⠂⡠⠂⠐⠐⠀⠀⠄⠀⠀⠀⠄⠀⠀⠀⠀⡀⠀⠀⠁⠀⠀⠀⠀⠀⠀⠠⠀⠂⠀⠠⠀⠀⠀⠆⠄⠈⠀⢀⠀⠀⠀⠐⠄⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⢠⠀⢀⠀⠀⠀⢀⠀⠄⠀⠀⠀⠂⠀⠀⠀⠇⠀⢀⠀⠀⠁⠁⠀⠄⠄⠀⠀⠁⠁⠀⡈⠀⠁⠄⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀│
      │⠀⠀⠀⠀⠀⠀⠀⠀⠀⡐⢀⠀⠀⠀⠐⠂⠁⢀⠠⠀⠀⠀⠀⠀⠁⠀⠀⠀⠠⠀⠀⠀⠀⢀⠀⢄⠀⠀⠀⠐⠀⠀⠈⠀⠀⠐⠀⠄⠂⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠠⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠠⠈⠀⠀⠀⠀⢀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
      │⠐⠀⠉⠐⠀⠀⠁⠀⠀⠀⠀⠀⠀⠌⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠌⠀⠀⠄⠀⠀⡀⠀⠀⠀⠀⠀⢀⠈⠄⠀⡀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀│
      │⠂⠀⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⡨⠀⠠⠀⠀⠀⠂⠀⠀⠀⠀⠀⠀⠈⡠⠀⡀⠀⠀⠁⠀⠁⠠⠀⠀⠀⠀⠀⠀⠀⠄⡀⠀⠄⡠⠀⠀⠀⠀│
      │⠐⠀⠀⠀⠀⠀⠀⠃⠀⠀⠀⠀⠁⠀⡀⠁⠀⢂⠂⠀⠀⠀⠈⠁⠀⠀⠁⠁⠄⠀⠀⠀⠀⠀⡀⠀⠀⠀⠂⠀⠀⠀⠀⠐⠀⠂⠀⠂⠁⠀⠆⠈⠀⠀⠄⠄⠀⠀⠀⠀│
      │⠀⠀⠀⠠⠀⠁⠀⠀⠀⠀⠂⠀⠠⠀⠀⠌⠀⠀⠂⠄⠀⡀⠂⠀⠀⡀⠠⠁⠄⠀⠀⠀⠀⠀⠄⠈⠀⢂⠐⠀⢀⠀⠠⠀⠄⠁⡀⠀⠀⠐⠀⠀⡐⠀⠀⠀⠁⠀⠀⠔│
      │⠀⠀⠀⠀⢀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠁⠃⡀⡁⠀⠀⠀⢀⠀⠀⠂⠀⠀⠠⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠬⠀⡀│
      │⠐⠀⠀⠀⠀⠀⠀⠀⠢⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⡀⠀⠀⠀⠀⠀⠈⠁⠁⠀⠀⠄⠀⠀⠐⠀⠄⠀⠠⠀⠁⡀⠁⠀⠀⠠⠡⠀⠀⠀⠀⠀⠀⠀⠐⠄⠀⠀⠀⠀⡀│
   50 │⠀⠀⠀⠀⠁⠀⠀⠠⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠂⠀⠀⠀⠀⠀⠀⠐⠀⠈⠂⠀⠀⠀⠀⠀⠀⠀⠀⠠⠈⠀⠀⠀⠀⠀⠀⠂⠀⠄⠀⠀⠀⠂⠁⠀⠀⠀⠂⠀│
      └────────────────────────────────────────────────────────────┘
      1                                                          120
                               nz = 275
```

# See also

[`Plot`](@ref), [`scatterplot`](@ref),
[`BrailleCanvas`](@ref), [`BlockCanvas`](@ref),
[`AsciiCanvas`](@ref), [`DotCanvas`](@ref)
"""
function spy(A::AbstractMatrix; kw...)
    rows, cols, vals = _strict_non_zeros(_findnz(A)...)
    if get(kw, :show_zeros, false)
        I = CartesianIndex.(zip(rows, cols))  # non zeros
        mask = trues(size(A))
        mask[I] .= false
        Z = CartesianIndices(axes(A))[mask]  # zeros
        rows, cols = getindex.(Z, 1), getindex.(Z, 2)
        vals = zeros(eltype(vals), length(rows))
    end
    spy(size(A)..., rows, cols, vals; kw...)
end

function _strict_non_zeros(rows, cols, vals)
    # findnz(A) returns stored zeros, ignore those
    I = findall(!iszero, vals)
    rows[I], cols[I], vals[I]
end

function spy(
    nrow::Int,
    ncol::Int,
    rows::AbstractArray{<:Integer},
    cols::AbstractArray{<:Integer},
    vals::AbstractArray;
    maxwidth::Int = 0,
    maxheight::Int = 0,
    title = "Sparsity Pattern",
    out_stream::Union{Nothing,IO} = nothing,
    width::Int = 0,
    height::Int = 0,
    margin::Int = KEYWORDS.margin,
    padding::Int = KEYWORDS.padding,
    color::UserColorType = KEYWORDS.color,
    canvas::Type{T} = BrailleCanvas,
    show_zeros::Bool = false,
    fix_ar::Bool = false,
    kw...,
) where {T<:Canvas}
    if color == :automatic
        Base.depwarn(
            "`color = :automatic` is deprecated, use `color = :auto` instead",
            :spy,
        )
        color = :auto
    end
    width, height = get_canvas_dimensions_for_matrix(
        canvas,
        nrow,
        ncol,
        maxwidth,
        maxheight,
        width,
        height,
        margin,
        padding,
        out_stream,
        fix_ar;
        extra_rows = 9,
        extra_cols = 6,
    )
    can = T(width, height, width = 1.0 + ncol, height = 1.0 + nrow)
    plot = Plot(can; title = title, margin = margin, padding = padding, kw...)

    if color != :auto
        points!(plot, cols, nrow + 1 .- rows, color)
        label!(plot, :r, 1, show_zeros ? "= 0" : "≠ 0", color)
    else
        if show_zeros
            points!(plot, cols, nrow + 1 .- rows, :green)
            label!(plot, :r, 1, "= 0", :green)
        else
            pos_idx = vals .> 0
            neg_idx = (!).(pos_idx)
            pos_cols = cols[pos_idx]
            pos_rows = rows[pos_idx]
            neg_cols = cols[neg_idx]
            neg_rows = rows[neg_idx]
            points!(plot, pos_cols, nrow + 1 .- pos_rows, :red)
            points!(plot, neg_cols, nrow + 1 .- neg_rows, :blue)
            label!(plot, :r, 1, "> 0", :red)
            label!(plot, :r, 2, "< 0", :blue)
        end
    end
    label!(plot, :l, 1, "1", :light_black)
    label!(plot, :l, nrows(plot.graphics), string(nrow), :light_black)
    label!(plot, :bl, "1", :light_black)
    label!(plot, :br, string(ncol), :light_black)
    haskey(kw, :xlabel) ||
        xlabel!(plot, string(length(vals), show_zeros ? " zeros" : " nonzeros"))
    plot
end

function _findnz(A::AbstractMatrix)
    I = findall(!iszero, A)
    getindex.(I, 1), getindex.(I, 2), A[I]
end

_findnz(A::AbstractSparseMatrix) = findnz(A)

ext: jl


text:
using Vectorize
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

## Test Apple Accelerate
@static if is_apple()
    include("AccelerateTests.jl")
else
    println("Accelerate not present. Aborting Accelerate tests")
end

## Test Yeppp
include("YepppTests.jl")

## Check for presence of VML and test VML
if Libdl.find_library(["libmkl_vml_avx"], ["/opt/intel/mkl/lib"]) != ""
    include("VMLTests.jl")
end

## Running tests over benchmarked Vectorized functions
include("VectorizeTests.jl")


text:
#=
property.jl
=#

using FillArrays: Zeros
#using ImageGeoms #: [many]

"""
    `ImageGeom` help

## Properties

Access by `ig.key`, mostly for 2D (and 3D) geometries:

`fovs` `(|dx|*nx, |dy|*ny, ...)`
`np`   `sum(mask)` = # pixels to be estimated

`dim`  `(nx, ny, [nz])`

`wx`   `(nx - 1/2) * dx + offset_x`
`wy`   `(ny - 1/2) * dy + offset_y`
`wz`   `(nz - 1/2) * dz + offset_z`

`x` 1D x coordinates of each pixel
`y` 1D y coordinates of each pixel
`z` 1D z coordinates of each pixel

`xg` x coordinates of each pixel as a grid (2D or 3D)
`yg` y coordinates of each pixel as a grid (2D or 3D)
`zg` z coordinates of each pixel as a grid (2D or 3D)

`u` 1D frequency domain coordinates of each pixel (kx)
`v` 1D frequency domain coordinates of each pixel (ky)
`w` 1D frequency domain coordinates of each pixel (kz)

`ug` ≥1D grid of kx frequency domain coordinates
`vg` ≥2D grid of ky frequency domain coordinates
`wg` ≥3D grid of kz frequency domain coordinates

`ones`  `FillArray.Ones(dim)`
`zeros` `FillArray.Zeros(dim)`

`mask_outline` binary image showing 2D outline of `mask` (todo)

Derived values for 3D images:

`mask_or` `[nx ny]` logical 'or' of `mask`

## Methods:

* `embed(x)`   turn short column(s) into array(s)
* `maskit(x)`  opposite of embed
* `shape(x)`   reshape long column x to `[nx ny [nz]]` array
* `unitv(...)`
  * `j=j | i=(ix,iy) | c=(cx cy)`
  * `j`: single index from 1 to `length(z)`
  * `i`: (ix,iy[,iz]) index from 1 to nx,ny
  * `c`: (cx,cy[,cz]) index from +/- n/2 center at floor(n/2)+1

* `circ(rx=,ry=,cx=,cy=)` circle of given radius and center (cylinder in 3D)
* `plot(jim)` plot the image geometry using the `MIRTjim.jim` function

## Methods that return a new `ImageGeom:`

* `down(down::Int)` down-sample geometry by given factor
* `over(over::Int)` over-sample geometry by given factor
* `expand_nz(nz_pad)` expand image geometry in z by `nz_pad` on both ends
"""
ImageGeomHelp

"""
    out = image_geom_unitv(dims::Dims ; j=?, i=?, c=?, T=?)

Return a standard unit vector (Kronecker impulse).

# options (use at one of these):
- `j` single index from 1 to length(z)
- `i` (ix, iy [,iz]) index from 1 to nx,ny
- `c` (cx, cy [,cz]) index from +/- n/2 center at floor(n/2)+1
- `T::DataType` default `Bool`

Default with no arguments gives unit vector at center `c=(0, 0 [,0])`
"""
function image_geom_unitv(
    dims::Dims{D} ;
    T::DataType = Bool,
    j::Int = 0,
    i::NTuple{D,Int} = ntuple(i -> 0, D),
    c::NTuple{D,Int} = ntuple(i -> 0, D),
)::Array{T, D} where {D}

    out = zeros(T, dims)
@show T, typeof(out)

#=
    (j < 0 || j > length(out)) && throw("bad j $j")
    if 1 <= j <= length(out)
        any(!=(0), i) && throw("i $i")
        any(!=(0), c) && throw("c $c")
        out[j] = one(T)
    elseif all(1 .<= i .<= size(out))
        any(!=(0), c) && throw("c $c")
        out[i...] = one(T)
    else
        any(!=(0), i) && throw("i $i")
        tmp = c .+ (size(out) .÷ 2) .+ 1
        out[tmp...] = one(T)
    end
=#

    return out::Array{T, D}
end

# spatial axes grids ala ndgrid
xg(ig::ImageGeom) = grids(ig)[1]
yg(ig::ImageGeom) = grids(ig)[2]
zg(ig::ImageGeom) = grids(ig)[3]

# spatial frequency axes grids ala ndgrid
ug(ig::ImageGeom) = gridf(ig)[1]
vg(ig::ImageGeom) = gridf(ig)[2]
wg(ig::ImageGeom) = gridf(ig)[3]

function help(ig::ImageGeom; io::IO = isinteractive() ? stdout : devnull)
    print(io, "propertynames:\n")
    println(io, propertynames(ig))
    @doc ImageGeomHelp
end

# Extended properties

image_geom_fun0 = Dict([
    (:help, ig -> help(ig)),

    (:ndim, ig -> length(ig.dims)),
    (:nx, ig -> ig.dims[1]),
    (:ny, ig -> ig.ndim ≥ 2 ? ig.dims[2] : 0),
    (:nz, ig -> ig.ndim ≥ 3 ? ig.dims[3] : 0),
    (:dim, ig -> ig.dims),
    (:fovs, ig -> abs.(ig.deltas) .* ig.dims),

    (:zeros, ig -> zeros(ig)),
    (:ones, ig -> ones(ig)),

    (:dx, ig -> ig.deltas[1]),
    (:dy, ig -> ig.ndim ≥ 2 ? ig.deltas[2] : zero(ig)),
    (:dz, ig -> ig.ndim ≥ 3 ? ig.deltas[3] : zero(ig)),

    (:offset_x, ig -> ig.offsets[1]),
    (:offset_y, ig -> ig.ndim ≥ 2 ? ig.offsets[2] : Float32(0)),
    (:offset_z, ig -> ig.ndim ≥ 3 ? ig.offsets[3] : Float32(0)),

    (:wx, ig -> (ig.nx - 1)/2 + ig.offset_x),
    (:wy, ig -> (ig.ny - 1)/2 + ig.offset_y),
    (:wz, ig -> (ig.nz - 1)/2 + ig.offset_z),

    # spatial axes and grids
    (:x, ig -> ((0:ig.nx-1) .- ig.wx) * ig.dx),
    (:y, ig -> ((0:ig.ny-1) .- ig.wy) * ig.dy),
    (:z, ig -> ((0:ig.nz-1) .- ig.wz) * ig.dz),

    (:xg, ig -> xg(ig)),
    (:yg, ig -> yg(ig)),
    (:zg, ig -> zg(ig)),

    # DFT frequency axes and grids
    (:u, ig -> axisf(ig, 1)),
    (:v, ig -> axisf(ig, 2)),
    (:w, ig -> axisf(ig, 3)),

    (:ug, ig -> ug(ig)),
    (:vg, ig -> vg(ig)),
    (:wg, ig -> wg(ig)),

    (:np, ig -> sum(ig.mask)),
    (:mask_or, ig -> mask_or(ig.mask)),
#   (:mask_outline, ig -> mask_outline(ig.mask)),

    # simple functions

    (:plot, ig -> ((how::Function ; kwargs...) -> plot(ig, how ; kwargs...))),
    (:embed, ig -> (x::AbstractArray -> embed(x, ig.mask))),
    (:maskit, ig -> (x::AbstractArray -> maskit(x, ig.mask))),
    (:shape, ig -> (x::AbstractArray -> reshape(x, ig.dims))),
    (:unitv, ig -> (( ; kwargs...) -> image_geom_unitv(ig.dim ; kwargs...))),
    (:circ, ig -> (( ; kwargs...) -> circle(ig ; kwargs...))),

    # functions that return new geometry

    (:down, ig -> (down::Int -> downsample(ig, down))),
    (:over, ig -> (over::Int -> oversample(ig, over))),
    (:expand_nz, ig -> (nz_pad::Int -> expand_nz(ig, nz_pad))),

    ]
)

# Tricky overloading here!

Base.getproperty(ig::ImageGeom, s::Symbol) =
        haskey(image_geom_fun0, s) ? image_geom_fun0[s](ig) : getfield(ig, s)

Base.propertynames(ig::ImageGeom) =
    (fieldnames(typeof(ig))..., keys(image_geom_fun0)...)


text:
"""
Read a tea header from the given IO stream.

`io` will be at the start of the item area following this function call, which could be
the end of the file.
"""
Base.read(io::IO, ::Type{TeaFileMetadata}) = _read_tea(io, TeaFileMetadata)

function _read_tea(io::IO, ::Type{TeaFileMetadata})::TeaFileMetadata
    _read_magic(io)

    item_start = read(io, Int64)
    item_end = read(io, Int64)
    n_sections = read(io, Int64)

    sections = AbstractSection[]
    for i_section in 1:n_sections
        push!(sections, _read_section(io))
    end

    pos = position(io)
    if pos > item_start
        error("Invalid tea-file, item_start=$item_start is within header.")
    elseif pos < item_start
        # Gobble up padding until the start of the item area.
        seek(io, item_start)
    end

    return TeaFileMetadata(item_start, item_end, sections)
end

"""Consume the magic value if present, or throw if not."""
function _read_magic(io::IO)
    marker = read(io, Int64)
    if marker != MAGIC_VALUE
        if marker == bswap(MAGIC_VALUE)
            throw(ArgumentError(
                "This is a tea-file written on a system with different endianness."
            ))
        else
            throw(ArgumentError(
                "Expected file to start with $MAGIC_VALUE, but got $marker."
            ))
        end
    end
end

"""Read a section of the correct type, and return an offset to the next section."""
function _read_section(io::IO)::AbstractSection
    id = read(io, Int32)
    next_section_offset = read(io, Int32)
    pos_before = position(io)
    section = _read_tea(io, section_type(id))
    pos_after = position(io)
    # Advance the cursor if required -- this will skip any padding that might exist between
    # this and the next section.
    skip(io, next_section_offset - (pos_after - pos_before))
    return section
end

"""Read a value previously written with _write_tea."""
_read_tea(io::IO, ::Type{T}) where T <: Real = read(io, T)
_read_tea(io::IO, ::Type{Base.UUID}) = Base.UUID(read(io, UInt128))

function _read_tea(io::IO, ::Type{String})::String
    length = read(io, Int32)
    return String(read(io, length))
end

function _read_tea(io::IO, ::Type{Vector{T}})::Vector{T} where T
    length = read(io, Int32)
    return [_read_tea(io, T) for _ in 1:length]
end

function _read_tea(io::IO, ::Type{Field})
    type_id = read(io, Int32)
    offset = read(io, Int32)
    name = _read_tea(io, String)
    return Field(type_id, offset, name)
end

function _read_tea(io::IO, ::Type{NameValue{T}})::NameValue{T} where T
    name_value = _read_tea(io, NameValue)
    type = typeof(name_value.value)
    if type != T
        error("$T does not match the type read, which is $type.")
    end
    return name_value
end

function _read_tea(io::IO, ::Type{NameValue})::NameValue
    name = _read_tea(io, String)
    kind = read(io, Int32)
    T = name_value_type(kind)
    value = _read_tea(io, T)
    return NameValue(name, value)
end

function _read_tea(io::IO, ::Type{ItemSection})::ItemSection
    item_size = read(io, Int32)
    item_name = _read_tea(io, String)
    fields = _read_tea(io, Vector{Field})
    return ItemSection(item_size, item_name, fields)
end

function _read_tea(io::IO, ::Type{TimeSection})::TimeSection
    epoch = read(io, Int64)
    ticks_per_day = read(io, Int64)
    time_field_offsets = _read_tea(io, Vector{Int32})
    return TimeSection(epoch, ticks_per_day, time_field_offsets)
end

function _read_tea(io::IO, ::Type{ContentDescriptionSection})::ContentDescriptionSection
    return ContentDescriptionSection(_read_tea(io, String))
end

function _read_tea(io::IO, ::Type{NameValueSection})::NameValueSection
    return NameValueSection(_read_tea(io, Vector{NameValue}))
end


text:
using AdvancedMH
using Distributions
using StructArrays
using MCMCChains
using LinearAlgebra

using Random
using Test
using DiffResults
using ForwardDiff

include("util.jl")

@testset "AdvancedMH" begin
    # Set a seed
    Random.seed!(1234)

    # Generate a set of data from the posterior we want to estimate.
    data = rand(Normal(0, 1), 300)

    # Define the components of a basic model.
    insupport(θ) = θ[2] >= 0
    dist(θ) = Normal(θ[1], θ[2])
    density(θ) = insupport(θ) ? sum(logpdf.(dist(θ), data)) : -Inf

    # Construct a DensityModel.
    model = DensityModel(density)

    @testset "StaticMH" begin
        # Set up our sampler with initial parameters.
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])
        spl2 = StaticMH(MvNormal([0.0, 0.0], 1))

        # Sample from the posterior.
        chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"])
        chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"])

        # chn_mean ≈ dist_mean atol=atol_v
        @test mean(chain1.μ) ≈ 0.0 atol=0.1
        @test mean(chain1.σ) ≈ 1.0 atol=0.1
        @test mean(chain2.μ) ≈ 0.0 atol=0.1
        @test mean(chain2.σ) ≈ 1.0 atol=0.1
    end

    @testset "RandomWalk" begin
        # Set up our sampler with initial parameters.
        spl1 = RWMH([Normal(0,1), Normal(0, 1)])
        spl2 = RWMH(MvNormal([0.0, 0.0], 1))

        # Sample from the posterior.
        chain1 = sample(model, spl1, 100000; chain_type=StructArray, param_names=["μ", "σ"])
        chain2 = sample(model, spl2, 100000; chain_type=StructArray, param_names=["μ", "σ"])

        # chn_mean ≈ dist_mean atol=atol_v
        @test mean(chain1.μ) ≈ 0.0 atol=0.1
        @test mean(chain1.σ) ≈ 1.0 atol=0.1
        @test mean(chain2.μ) ≈ 0.0 atol=0.1
        @test mean(chain2.σ) ≈ 1.0 atol=0.1
    end

    @testset "AdaptiveMH" begin
        # Set up our sampler with initial parameters.
        doubledensity(θ::Vector{Float64}) = density(θ[1:2]) + density(θ[3:4])
        doubledensity(θ::Vector{Vector{Float64}}) = density(θ[1]) + density(θ[2])
        doublemodel = DensityModel(doubledensity)

        spl1 = MetropolisHastings(AMProposal(diagm(0 => [0.01, 0.01, 0.01, 0.01])))
        spl2 = MetropolisHastings([AMProposal(diagm(0 => [0.01, 0.01])),
                                   AMProposal(diagm(0 => [0.01, 0.01]))])

        # Sample from the posterior.
        chain1 = sample(doublemodel, spl1, 100000; chain_type=StructArray,
                        param_names=["μ1", "σ1", "μ2", "σ2"])
        chain2 = sample(doublemodel, spl2, 100000; chain_type=StructArray,
                        param_names=["p1", "p2"])

        # chn_mean ≈ dist_mean atol=atol_v
        @test mean(chain1.μ1) ≈ 0.0 atol=0.1
        @test mean(chain1.σ1) ≈ 1.0 atol=0.1
        @test mean(chain1.μ2) ≈ 0.0 atol=0.1
        @test mean(chain1.σ2) ≈ 1.0 atol=0.1

        @test mean(chain2.p1) ≈ [ 0.0 1.0 ]' atol=0.1
        @test mean(chain2.p2) ≈ [ 0.0 1.0 ]' atol=0.1
    end

    @testset "parallel sampling" begin
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])

        chain1 = sample(model, spl1, MCMCDistributed(), 10000, 4;
                        param_names=["μ", "σ"], chain_type=Chains)
        @test mean(chain1["μ"]) ≈ 0.0 atol=0.1
        @test mean(chain1["σ"]) ≈ 1.0 atol=0.1

        if VERSION >= v"1.3"
            chain2 = sample(model, spl1, MCMCThreads(), 10000, 4;
                            param_names=["μ", "σ"], chain_type=Chains)
            @test mean(chain2["μ"]) ≈ 0.0 atol=0.1
            @test mean(chain2["σ"]) ≈ 1.0 atol=0.1
        end
    end

    @testset "MCMCChains" begin
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])
        spl2 = MetropolisHastings((μ = StaticProposal(Normal(0,1)), σ = StaticProposal(Normal(0, 1))))

        chain1 = sample(model, spl1, 10_000; param_names=["μ", "σ"], chain_type=Chains)
        chain2 = sample(model, spl2, 10_000; chain_type=Chains)

        @test mean(chain1["μ"]) ≈ 0.0 atol=0.1
        @test mean(chain1["σ"]) ≈ 1.0 atol=0.1

        @test mean(chain2["μ"]) ≈ 0.0 atol=0.1
        @test mean(chain2["σ"]) ≈ 1.0 atol=0.1
    end

    @testset "Proposal styles" begin
        m1 = DensityModel(x -> logpdf(Normal(x,1), 1.0))
        m2 = DensityModel(x -> logpdf(Normal(x[1], x[2]), 1.0))
        m3 = DensityModel(x -> logpdf(Normal(x.a, x.b), 1.0))
        m4 = DensityModel(x -> logpdf(Normal(x,1), 1.0))

        p1 = StaticProposal(Normal(0,1))
        p2 = StaticProposal([Normal(0,1), InverseGamma(2,3)])
        p3 = (a=StaticProposal(Normal(0,1)), b=StaticProposal(InverseGamma(2,3)))
        p4 = StaticProposal((x=1.0) -> Normal(x, 1))

        c1 = sample(m1, MetropolisHastings(p1), 100; chain_type=Vector{NamedTuple})
        c2 = sample(m2, MetropolisHastings(p2), 100; chain_type=Vector{NamedTuple})
        c3 = sample(m3, MetropolisHastings(p3), 100; chain_type=Vector{NamedTuple})
        c4 = sample(m4, MetropolisHastings(p4), 100; chain_type=Vector{NamedTuple})

        @test keys(c1[1]) == (:param_1, :lp)
        @test keys(c2[1]) == (:param_1, :param_2, :lp)
        @test keys(c3[1]) == (:a, :b, :lp)
        @test keys(c4[1]) == (:param_1, :lp)
    end
    @testset "Initial parameters" begin
        # Set up our sampler with initial parameters.
        spl1 = StaticMH([Normal(0,1), Normal(0, 1)])

        val = [0.4, 1.2]

        # Sample from the posterior.
        chain1 = sample(model, spl1, 10, init_params = val)

        @test chain1[1].params == val
    end

    @testset "symmetric proposals" begin
        # True distributions
        d1 = Normal(5, .7)

        # Model definition.
        m1 = DensityModel(x -> logpdf(d1, x))

        # Custom normal distribution without `logpdf` defined errors since the
        # acceptance probability cannot be computed
        p1 = RandomWalkProposal(CustomNormal())
        @test p1 isa RandomWalkProposal{false}
        @test_throws MethodError AdvancedMH.logratio_proposal_density(p1, randn(), randn())
        @test_throws MethodError sample(m1, MetropolisHastings(p1), 10)

        p1 = StaticProposal(x -> CustomNormal(x))
        @test p1 isa StaticProposal{false}
        @test_throws MethodError AdvancedMH.logratio_proposal_density(p1, randn(), randn())
        @test_throws MethodError sample(m1, MetropolisHastings(p1), 10)

        # If the proposal is declared to be symmetric, the log ratio of the proposal
        # density is not evaluated.
        p2 = SymmetricRandomWalkProposal(CustomNormal())
        @test p2 isa RandomWalkProposal{true}
        p2 = SymmetricStaticProposal(x -> CustomNormal(x))
        @test p2 isa StaticProposal{true}

        for p2 in (
            SymmetricRandomWalkProposal(CustomNormal()),
            SymmetricStaticProposal((x=0) -> CustomNormal(x)),
        )
            @test iszero(AdvancedMH.logratio_proposal_density(p2, randn(), randn()))
            @test iszero(AdvancedMH.logratio_proposal_density([p2], randn(1), randn(1)))
            @test iszero(AdvancedMH.logratio_proposal_density(
                (p2,), (randn(),), (randn(),)
            ))
            @test iszero(AdvancedMH.logratio_proposal_density(
                (; x=p2), (; x=randn()), (; x=randn())
            ))
            chain1 = sample(
                m1, MetropolisHastings(p2), 100000;
                chain_type=StructArray, param_names=["x"]
            )
            @test mean(chain1.x) ≈ mean(d1) atol=0.05
            @test std(chain1.x) ≈ std(d1) atol=0.05
        end

        # type inference checks (arrays of proposals are not guaranteed to be type stable)
        proposals = (
            StaticProposal(Normal()),
            StaticProposal(x -> Normal(x, 1)),
            StaticProposal{true}(Cauchy()),
            StaticProposal{true}(x -> Cauchy(x, 2)),
            RandomWalkProposal(Laplace()),
            RandomWalkProposal(x -> Laplace(x, 1)),
            RandomWalkProposal{true}(TDist(1)),
            RandomWalkProposal{true}(x -> TDist(1)),
        )
        states = randn(2)
        candidates = randn(2)
        for (p1, p2) in Iterators.product(proposals, proposals)
            val = AdvancedMH.logratio_proposal_density(p1, states[1], candidates[1]) +
                AdvancedMH.logratio_proposal_density(p2, states[2], candidates[2])
            @test AdvancedMH.logratio_proposal_density([p1, p2], states, candidates) ≈ val
            @test @inferred AdvancedMH.logratio_proposal_density(
                (p1, p2), (states[1], states[2]), (candidates[1], candidates[2])
            ) ≈ val
            @test @inferred AdvancedMH.logratio_proposal_density(
                (x=p1, y=p2), (y=states[2], x=states[1]), (x=candidates[1], y=candidates[2])
            ) ≈ val
        end
    end

    @testset "MALA" begin

        # Set up the sampler.
        sigma = 1e-1
        spl1 = MALA(x -> MvNormal((sigma^2 / 2) .* x, sigma))

        # Sample from the posterior with initial parameters.
        chain1 = sample(model, spl1, 100000; init_params=ones(2), chain_type=StructArray, param_names=["μ", "σ"])

        @test mean(chain1.μ) ≈ 0.0 atol=0.1
        @test mean(chain1.σ) ≈ 1.0 atol=0.1
    end

    @testset "EMCEE" begin include("emcee.jl") end
end


text:
struct KnowledgeNetworkforBiocomplexity <: KNB
end

function pub_date(repo::KnowledgeNetworkforBiocomplexity, mainpage)
    "No dates released"
end


text:
using POMDPs
using POMDPToolbox
using POMDPModels
using Base.Test

include("test_random_solver.jl")
include("test_stochastic_policy.jl")
include("test_rollout.jl")
include("test_history_recorder.jl")
include("test_sim.jl")
include("test_stepthrough.jl")
include("test_parallel.jl")
include("test_belief.jl")
include("test_particle.jl")
include("test_solver_test.jl")
include("test_model.jl")
include("test_generative_belief_mdp.jl")
include("test_alpha_policy.jl")
include("test_vector_policy.jl")
include("test_value_policy.jl")
include("test_function_policy.jl")
include("test_utility_wrapper.jl")
include("test_implementations.jl")
include("test_distributions_jl.jl")
include("test_weighted_iteration.jl")
include("test_sparse_cat.jl")
include("test_bool.jl")
include("test_info.jl")
include("test_k_previous_observations_belief.jl")
include("test_fully_observable_pomdp.jl")
include("test_underlying_mdp.jl")
include("test_obs_weight.jl")


text:
using Pkg; Pkg.activate(@__DIR__); Pkg.instantiate()

using Trixi

###############################################################################
# Implementation of the flux of Sjögreen and Yee (2019) for the compressible
# Euler equations on the 2D `TreeMesh` of Trixi.jl.
# https://doi.org/10.1007/s10915-019-01013-1

struct HartensCompressibleEulerEquations2D{RealT<:Real} <: Trixi.AbstractCompressibleEulerEquations{2, 4}
  equations::CompressibleEulerEquations2D{RealT}
  alpha::RealT # parameter of Harten's entropies
end

function HartensCompressibleEulerEquations2D(gamma::Real, alpha::Real)
  equations = CompressibleEulerEquations2D(gamma)
  HartensCompressibleEulerEquations2D(equations, alpha)
end

Trixi.varnames(f, equations::HartensCompressibleEulerEquations2D) = Trixi.varnames(f, equations.equations)

Trixi.flux(u, orientation_or_normal_direction, equations::HartensCompressibleEulerEquations2D) = flux(u, orientation_or_normal_direction, equations.equations)

Trixi.flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction, equations::HartensCompressibleEulerEquations2D) = flux_shima_etal(u_ll, u_rr, orientation_or_normal_direction, equations.equations)

Trixi.flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, equations::HartensCompressibleEulerEquations2D) = flux_ranocha(u_ll, u_rr, orientation_or_normal_direction, equations.equations)

Trixi.initial_condition_density_wave(x, t, equations::HartensCompressibleEulerEquations2D) = initial_condition_density_wave(x, t, equations.equations)

Trixi.initial_condition_weak_blast_wave(x, t, equations::HartensCompressibleEulerEquations2D) = initial_condition_weak_blast_wave(x, t, equations.equations)

Trixi.initial_condition_convergence_test(x, t, equations::HartensCompressibleEulerEquations2D) = initial_condition_convergence_test(x, t, equations.equations)

Trixi.source_terms_convergence_test(u, x, t, equations::HartensCompressibleEulerEquations2D) = source_terms_convergence_test(u, x, t, equations.equations)

Trixi.cons2prim(u, equations::HartensCompressibleEulerEquations2D) = cons2prim(u, equations.equations)

Trixi.max_abs_speeds(u, equations::HartensCompressibleEulerEquations2D) = Trixi.max_abs_speeds(u, equations.equations)

Trixi.max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations::HartensCompressibleEulerEquations2D) = Trixi.max_abs_speed_naive(u_ll, u_rr, orientation_or_normal_direction, equations.equations)

function Trixi.entropy(u, equations::HartensCompressibleEulerEquations2D)
  alpha = equations.alpha
  gamma = equations.equations.gamma
  rho, v1, v2, p = cons2prim(u, equations)

  return -(gamma + alpha) / (gamma - 1) * rho * (p / rho^gamma)^(1 / (alpha + gamma))
end

function Trixi.cons2entropy(u, equations::HartensCompressibleEulerEquations2D)
  alpha = equations.alpha
  gamma = equations.equations.gamma
  rho, v1, v2, p = cons2prim(u, equations)

  h = (p / rho^gamma)^(1 / (alpha + gamma))
  factor = rho / p * h
  w1 = -alpha / (gamma - 1) * h - 0.5 * factor * (v1^2 + v2^2)
  w2 = factor * v1
  w3 = factor * v2
  w5 = -factor

  return SVector(w1, w2, w3, w5)
end

function flux_potential(u, orientation::Integer, equations::HartensCompressibleEulerEquations2D)
  alpha = equations.alpha
  gamma = equations.equations.gamma
  rho, v1, v2, p = cons2prim(u, equations)

  h = (p / rho^gamma)^(1 / (alpha + gamma))

  if orientation == 1
    return rho * h * v1
  else # if orientation == 2
    return rho * h * v2
  end
end

function cons2param(u, equations::HartensCompressibleEulerEquations2D)
  alpha = equations.alpha
  gamma = equations.equations.gamma
  rho, v1, v2, p = cons2prim(u, equations)

  h = (p / rho^gamma)^(1 / (alpha + gamma))
  factor = rho / p * h
  z1 = factor
  z2 = v1
  z3 = v2
  z5 = p

  return SVector(z1, z2, z3, z5)
end

function exponential_average(a, b, beta_minus_1)
  beta = beta_minus_1 + 1

  r = a / b
  if abs(r - 1) < 1.0e-4
    return b^beta_minus_1 * (
      1 + (beta - 1) * (r - 1) / 2 + (beta - 1) * (beta - 2) * (r - 1)^2 / 6
    )
  else
    return (a^beta - b^beta) / (beta * (a - b))
  end
end

function flux_sjögreen_yee(u_ll, u_rr, orientation::Integer, equations::HartensCompressibleEulerEquations2D)
  alpha = equations.alpha
  gamma = equations.equations.gamma

  # Unpack left and right state
  z1_ll, z2_ll, z3_ll, z5_ll = cons2param(u_ll, equations)
  z1_rr, z2_rr, z3_rr, z5_rr = cons2param(u_rr, equations)

  if orientation == 1
    zhat_ll = z2_ll
    zhat_rr = z2_rr
  else # if orientation == 2
    zhat_ll = z3_ll
    zhat_rr = z3_rr
  end

  z1_zhat_mean = 0.5 * (z1_ll * zhat_ll + z1_rr * zhat_rr)
  z1_power_mean = 0.5 * (z1_ll^(-gamma / alpha) + z1_rr^(-gamma / alpha))
  z5_power_expmean = exponential_average(z5_ll, z5_rr, (1 - gamma - alpha) / alpha)
  h1 = z1_zhat_mean / (z1_power_mean * z5_power_expmean)

  if orientation == 1
    h2 = h1 * 0.5 * (z2_ll + z2_rr) + 0.5 * (z5_ll + z5_rr)
    h3 = h1 * 0.5 * (z3_ll + z3_rr)
  else # if orientation == 2
    h2 = h1 * 0.5 * (z2_ll + z2_rr)
    h3 = h1 * 0.5 * (z3_ll + z3_rr) + 0.5 * (z5_ll + z5_rr)
  end

  z5_power_mean = 0.5 * (z5_ll^(-(gamma - 1) / alpha) + z5_rr^(-(gamma - 1) / alpha))
  z1_power_expmean = exponential_average(z1_ll, z1_rr, -(gamma + alpha) / alpha)
  h5 = h1 * (
    gamma / (gamma - 1) * z5_power_mean * z1_power_expmean
    + 0.25 * (z2_ll + z2_rr)^2
    + 0.25 * (z3_ll + z3_rr)^2
    - 0.25 * (z2_ll^2 + z2_rr^2 + z3_ll^2 + z3_rr^2)
  )

  return SVector(h1, h2, h3, h5)
end


text:
# RCP3PD, including land use, GtC
# 270 (natural) + 352.34 (anthro) MtCH4
# 10.7 (natural) + 12.66*28/44 (anthro) MtN
# 		      breakyear1 breakyear2 growth1 growth2 growth3
const scendata = [
"BAU" 		:CO2	2050	2150	0.015	0.0018	0.0018
"BAU" 		:CH4	2050	2150	0.008	0		0
"BAU" 		:N2O	2050	2150	0.005	0		0
"3DEG" 		:CO2	2050	2150	-0.0067	-0.0067	-0.0067
"3DEG" 		:CH4	2050	2150	0		0		0
"3DEG" 		:N2O	2050	2150	0		0		0
"2DEG" 		:CO2	2050	2150	-0.0285	-0.0075	-0.0075
"2DEG" 		:CH4	2050	2150	0		0		0
"2DEG" 		:N2O	2050	2150	0		0		0
"1DEG" 		:CO2	2050	2150	-0.0617	-0.0022	-0.0022
"1DEG" 		:CH4	2050	2150	-0.01	0		0
"1DEG" 		:N2O	2050	2150	-0.01	0		0
"15DEG" 	:CO2	2050	2150	-0.055	-0.055	-0.055
"15DEG" 	:CH4	2050	2150	-0.0014	-0.0014	-0.0014
"15DEG" 	:N2O	2050	2150	-0.0014	-0.0014	-0.0014
"1PCT" 		:CO2	2050	2150	-0.01	-0.01	-0.01
"1PCT" 		:CH4	2050	2150	0		0		0
"1PCT" 		:N2O	2050	2150	0		0		0
"1PCTD" 	:CO2	2050	2150	0		-0.025	-0.025
"1PCTD" 	:CH4	2050	2150	0		-0.0005	-0.0005
"1PCTD" 	:N2O	2050	2150	0		-0.0005	-0.0005
"3GTC" 		:CO2	2050	2150	0.005	0		0
"3GTC" 		:CH4	2050	2150	0		0		0
"3GTC" 		:N2O	2050	2150	0		0		0
"2C_3CS"	:CO2	2100	2150	-0.0285	-0.015	-0.015
"2C_3CS"	:CH4	2050	2150	0		0		0
"2C_3CS"	:N2O	2050	2150	0		0		0
"2C_45CS"	:CO2	2100	2150	-0.0285	-0.015	-0.015
"2C_45CS"	:CH4	2100	2150	-0.0006	-0.0010	-0.0010
"2C_45CS"	:N2O	2100	2150	-0.0006	-0.0010	-0.0010
"2C_15CS"	:CO2	2065	2170	0.01	-0.019	-0.006
"2C_15CS"	:CH4	2050	2150	0		0		0
"2C_15CS"	:N2O	2050	2150	0		0		0
]

function getscenario(scen)
	breakyear1 = Dict{Tuple{String,Symbol},Int64}(
		(scendata[r,1], scendata[r,2]) => scendata[r,3] for r=1:size(scendata,1))
	breakyear2 = Dict{Tuple{String,Symbol},Int64}(
		(scendata[r,1], scendata[r,2]) => scendata[r,4] for r=1:size(scendata,1))
	growth1 = Dict{Tuple{String,Symbol},Float64}(
		(scendata[r,1], scendata[r,2]) => scendata[r,5] for r=1:size(scendata,1))
	growth2 = Dict{Tuple{String,Symbol},Float64}(
		(scendata[r,1], scendata[r,2]) => scendata[r,6] for r=1:size(scendata,1))
	growth3 = Dict{Tuple{String,Symbol},Float64}(
		(scendata[r,1], scendata[r,2]) => scendata[r,7] for r=1:size(scendata,1))

	path = joinpath(dirname(@__FILE__), "..")
	rcpscen = scen[1:3] == "RCP" ? scen : "RCP3PD"
	emissionsRCP = readfixedwidth("$path/RCP/$(rcpscen)_EMISSIONS.DAT", header=true, skipstart=37)

	#emissions2010 = Dict(:CO2 => 9.8779, :CH4 => 622.34, :N2O => 18.76)
	annualemissions = Dict{Symbol, Vector{Float64}}()
	annualemissions[:CO2] = emissionsRCP[:FossilCO2] + emissionsRCP[:OtherCO2]
	annualemissions[:CH4] = emissionsRCP[:CH4] .+ 270	# natural background emissions: 270 MtCH4
	annualemissions[:N2O] = emissionsRCP[:N2O] .+ 10.7	# natural background emissions: 10.7 MtN (TAR p253)
	if scen[1:3] != "RCP"
		for g in GAS3, y in 2010:YEARS[end]
			annualemissions[g][iyear(y)] = annualemissions[g][iyear(y)-1] *
				if y <= breakyear1[scen,g]
					1 + growth1[scen,g]
				elseif y <= breakyear2[scen,g]
					1 + growth2[scen,g]
				else
					1 + growth3[scen,g]
				end
		end
	end
	return annualemissions
end

max_stars_repo_name: niclasmattsson/ClimateCalculator

max_issues_repo_name: niclasmattsson/ClimateCalculator

max_forks_repo_name: niclasmattsson/ClimateCalculator


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for the Blade type.
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
import LinearAlgebra
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- Constructor tests

@testset "Blade: inner constructor" begin
    #=
      Notes
      -----
      * Test value of constructed instance
    =#

    # --- Preparations

    # Select random volume for tests where `volume` != 0
    test_volume = rand() + 1  # add 1 to avoid 0
    test_volume = rand() > 0.5 ? test_volume : -test_volume

    # --- Default constructor with constraint enforcement
    #
    # Blade{T}(dim::Int, grade::Int, basis::Matrix{T}, volume::Real;
    #          atol::Real=blade_atol(T), enforce_constraints::Bool=true,
    #          copy_basis::Bool=true) where {T<:AbstractFloat}

    # Preparations
    test_dim = 3
    test_grade = 2
    test_volume = 5

    for precision_type in subtypes(AbstractFloat)
        # --- Preparations

        test_basis = Matrix{precision_type}([3 -3; 4 -4; 0 1])
        normalized_basis = copy(test_basis)
        for i in 1:size(normalized_basis, 2)
            normalized_basis[:, i] /= LinearAlgebra.norm(normalized_basis[:, i])
        end

        # --- Valid data fields

        # volume > atol
        # default values for atol, enforce_constraints, copy_basis
        B = Blade{precision_type}(test_dim, test_grade, normalized_basis,
                                  test_volume)

        @test B.dim == test_dim
        @test B.grade == test_grade
        @test B.basis == normalized_basis
        @test B.basis !== normalized_basis
        for i in 1:size(B.basis, 2)
            @test LinearAlgebra.norm(B.basis[:, i]) ≈ 1
        end
        @test B.volume == test_volume

        # volume < atol
        # default values for enforce_constraints, copy_basis
        B = Blade{precision_type}(test_dim, test_grade, normalized_basis,
                                  test_volume, atol=abs(test_volume) + 1)
        @test iszero(B)

        # dim == 0
        B = Blade{precision_type}(0, test_grade,
                                  Matrix{precision_type}(rand(0, test_grade)),
                                  test_volume)
        @test iszero(B)

        # grade == 0
        B = Blade{precision_type}(test_dim, 0,
                                  Matrix{precision_type}(rand(test_dim, 0)),
                                  test_volume)
        @test iszero(B)

        # dim == grade, det(basis) > 0
        # default values for atol, enforce_constraints, copy_basis
        pseudoscalar_basis = Matrix{precision_type}([2 1; 4 3])
        for i in 1:size(pseudoscalar_basis, 2)
            pseudoscalar_basis[:, i] /=
                LinearAlgebra.norm(pseudoscalar_basis[:, i])
        end
        B = Blade{precision_type}(2, 2, pseudoscalar_basis, test_volume)

        @test B isa Pseudoscalar{precision_type}
        @test B.dim == 2
        @test B.value == test_volume

        # dim == grade, det(basis) < 0
        # default values for atol, enforce_constraints, copy_basis
        pseudoscalar_basis = Matrix{precision_type}([1 2; 3 4])
        for i in 1:size(pseudoscalar_basis, 2)
            pseudoscalar_basis[:, i] /=
                LinearAlgebra.norm(pseudoscalar_basis[:, i])
        end
        B = Blade{precision_type}(2, 2, pseudoscalar_basis, test_volume)

        @test B isa Pseudoscalar{precision_type}
        @test B.dim == 2
        @test B.value == -test_volume

        # dim == grade, det(basis) == 0
        # default values for atol, enforce_constraints, copy_basis
        pseudoscalar_basis = Matrix{precision_type}([1 4; 2 8])
        for i in 1:size(pseudoscalar_basis, 2)
            pseudoscalar_basis[:, i] /=
                LinearAlgebra.norm(pseudoscalar_basis[:, i])
        end
        B = Blade{precision_type}(2, 2, pseudoscalar_basis, test_volume)
        @test iszero(B)

        # copy_basis = false
        B = Blade{precision_type}(test_dim, test_grade, normalized_basis,
                                  test_volume, copy_basis=false)

        @test B.basis === normalized_basis

        # --- Invalid data fields

        # ------ enforce_constraints = true

        # dim < 0
        @test_throws ArgumentError Blade{precision_type}(-5,
                                                          test_grade,
                                                          test_basis,
                                                          test_volume)

        # grade < 0
        @test_throws ArgumentError Blade{precision_type}(test_dim,
                                                          -1,
                                                          test_basis,
                                                          test_volume)

        # dim != size(basis, 1)
        @test_throws DimensionMismatch Blade{precision_type}(test_dim + 1,
                                                             test_grade,
                                                             test_basis,
                                                             test_volume)

        # grade != size(basis, 2)
        @test_throws DimensionMismatch Blade{precision_type}(test_dim,
                                                             test_grade + 1,
                                                             test_basis,
                                                             test_volume)

        # basis not normalized
        @test_throws ArgumentError Blade{precision_type}(test_dim,
                                                          test_grade,
                                                          test_basis,
                                                          test_volume)

        # ------ enforce_constraints = false

        # dim < 0
        B = Blade{precision_type}(-5, test_grade, test_basis, test_volume,
                                  enforce_constraints=false)
        @test B.dim == -5

        # dim != size(basis, 1)
        # grade != size(basis, 2)
        # basis not normalized
        B = Blade{precision_type}(test_dim + 1, test_grade + 1, test_basis,
                                  test_volume, enforce_constraints=false)
        @test B.dim == test_dim + 1
        @test B.grade == test_grade + 1
        @test B.basis == test_basis
        @test B.basis !== test_basis
    end
end

@testset "Blade: outer constructors - Blade{T}(vectors::Matrix)" begin
    #=
      Notes
      -----
      * Test value of constructed instance.

      * Test behavior of keyword arguments: `volume`, `atol`
    =#

    # --- Preparations

    vectors = Matrix([3 3; -4 -4; 0 1])
    test_dim = 3
    test_volume = 5

    # --- Tests

    for precision_type in subtypes(AbstractFloat)

        # --- vectors::Matrix{<:AbstractFloat}

        for value_type in subtypes(AbstractFloat)
            # Preparations
            converted_vectors = convert(Matrix{value_type}, vectors)

            # number of vectors <= dimension of column space, volume == nothing
            B = Blade{precision_type}(converted_vectors)
            @test B isa Blade{precision_type}
            @test B.dim == 3
            @test B.grade == 2
            @test size(B.basis) == (3, 2)
            for i in 1:size(B.basis, 2)
                @test LinearAlgebra.norm(B.basis[:, i]) ≈ 1
            end
            @test B.volume ≈ 5

            # number of vectors <= dimension of column space, volume != nothing
            B = Blade{precision_type}(converted_vectors, volume=test_volume)
            @test B isa Blade{precision_type}
            @test B.dim == 3
            @test B.grade == 2
            @test size(B.basis) == (3, 2)
            for i in 1:size(B.basis, 2)
                @test LinearAlgebra.norm(B.basis[:, i]) ≈ 1
            end

            F = LinearAlgebra.qr(converted_vectors)
            signed_norm = prod(LinearAlgebra.diag(F.R))
            @test B.volume == sign(signed_norm) * precision_type(test_volume)

            # number of vectors <= dimension of column space
            # volume < default atol
            B = Blade{precision_type}(converted_vectors,
                                      volume=blade_atol(precision_type) / 2)
            @test iszero(B)

            # norm(blade) < atol
            B = Blade{precision_type}(converted_vectors, atol=6)
            @test iszero(B)

            # one of the dimensions of basis is 0
            @test Blade{precision_type}(rand(0, test_dim)) ==
                zero(Blade{precision_type})
            @test Blade{precision_type}(rand(test_dim, 0)) ==
                zero(Blade{precision_type})

            # number of vectors > dimension of column space
            special_case_vectors = Matrix{value_type}([1 2 3; 4 5 6])
            B = Blade{precision_type}(special_case_vectors)
            @test iszero(B)

            # vectors are linearly dependent, volume == nothing
            special_case_vectors = Matrix{value_type}([1 2 1; 1 2 4; 1 2 9])
            B = Blade{precision_type}(special_case_vectors)
            @test iszero(B)

            # vectors are linearly dependent, volume != nothing
            special_case_vectors = Matrix{value_type}([1 2 1; 1 2 4; 1 2 9])
            B = Blade{precision_type}(special_case_vectors, volume=test_volume)
            @test iszero(B)

            # number of vectors == dimension of column space, volume == nothing
            # det(vectors) > 0
            special_case_vectors = Matrix{value_type}([3 5; 4 12])
            B = Blade{precision_type}(special_case_vectors)
            @test B isa Pseudoscalar{precision_type}
            @test B.dim == 2
            @test B.value ≈ LinearAlgebra.det(special_case_vectors)

            # number of vectors == dimension of column space, volume == nothing
            # det(vectors) < 0
            special_case_vectors = Matrix{value_type}([5 3; 12 4])
            B = Blade{precision_type}(special_case_vectors)
            @test B isa Pseudoscalar{precision_type}
            @test B.dim == 2
            @test B.value ≈ LinearAlgebra.det(special_case_vectors)

            # number of vectors == dimension of column space, volume == nothing
            # det(vectors) == 0
            special_case_vectors = Matrix{value_type}([5 1; 5 1])
            B = Blade{precision_type}(special_case_vectors)
            @test iszero(B)

            # number of vectors == dimension of column space, volume != nothing
            special_case_vectors = Matrix{value_type}([3 5; 4 12])
            B = Blade{precision_type}(special_case_vectors, volume=test_volume)
            @test B isa Pseudoscalar{precision_type}
            @test B.dim == 2
            sign_basis = sign(LinearAlgebra.det(convert(Matrix{precision_type},
                                                        special_case_vectors)))
            @test B.value ≈ sign_basis * precision_type(test_volume)
        end

        # --- vectors::Matrix{<:Integer}
        #
        # Notes
        # -----
        # * We only need to verify the correctness of the type of the result
        #   for a few permutations of the arguments because the computations
        #   for the data fields are independent of the type of the `vectors`
        #   argument. The correctness of the data field computations is
        #   verified in the unit tests for the cases when `vectors` has type
        #   Matrix{<:AbstractFloat}.

        # subtypes(Signed)
        for value_type in subtypes(Signed)
            # Preparations
            converted_vectors = convert(Matrix{value_type}, vectors)

            # number of vectors <= dimension of column space
            B = Blade{precision_type}(converted_vectors)
            @test B isa Blade{precision_type}

            # number of vectors <= dimension of column space
            # volume < default atol
            B = Blade{precision_type}(converted_vectors,
                                      volume=blade_atol(precision_type) / 2)
            @test iszero(B)

            # vectors are linearly dependent
            special_case_vectors = Matrix{value_type}([1 2 1; 1 2 4; 1 2 9])
            B = Blade{precision_type}(special_case_vectors)
            @test iszero(B)

            # number of vectors == dimension of column space
            special_case_vectors = Matrix{value_type}([3 5; 4 12])
            B = Blade{precision_type}(special_case_vectors)
            @test B isa Pseudoscalar{precision_type}
        end

        # subtypes(Unsigned)
        for value_type in subtypes(Unsigned)
            # Preparations
            converted_vectors = convert(Matrix{value_type}, abs.(vectors))

            # number of vectors <= dimension of column space
            B = Blade{precision_type}(converted_vectors)
            @test B isa Blade{precision_type}

            # number of vectors <= dimension of column space
            # volume < default atol
            B = Blade{precision_type}(converted_vectors,
                                      volume=blade_atol(precision_type) / 2)
            @test iszero(B)

            # vectors are linearly dependent
            special_case_vectors = Matrix{value_type}([1 2 1; 1 2 4; 1 2 9])
            B = Blade{precision_type}(special_case_vectors)
            @test iszero(B)

            # number of vectors == dimension of column space
            special_case_vectors = Matrix{value_type}([3 5; 4 12])
            B = Blade{precision_type}(special_case_vectors)
            @test B isa Pseudoscalar{precision_type}
        end

        # --- Bool

        # Preparations
        converted_vectors = Matrix{Bool}(vectors .!= 0)

        # number of vectors <= dimension of column space
        B = Blade{precision_type}(converted_vectors)
        @test B isa Blade{precision_type}

        # number of vectors <= dimension of column space
        # volume < default atol
        B = Blade{precision_type}(converted_vectors,
                                  volume=blade_atol(precision_type) / 2)
        @test iszero(B)

        # vectors are linearly dependent
        special_case_vectors = Matrix{Bool}([1 0 1; 0 1 1; 1 0 1])
        B = Blade{precision_type}(special_case_vectors)
        @test iszero(B)

        # number of vectors == dimension of column space
        special_case_vectors = Matrix{Bool}([1 1; 0 1])
        B = Blade{precision_type}(special_case_vectors)
        @test B isa Pseudoscalar{precision_type}
    end
end

@testset "Blade: outer constructors - Blade{T}(v::Vector)" begin
    #=
      Notes
      -----
      * Test value of constructed instance.

      * Test behavior of keyword arguments: `volume`, `atol`
    =#
    # --- Preparations

    vector = [3, 4, 12]
    test_volume = 5

    # --- Tests

    for precision_type in subtypes(AbstractFloat)
        # Preparations

        # --- vector::Vector{<:AbstractFloat}

        for value_type in subtypes(AbstractFloat)
            # Preparations
            converted_vector = convert(Vector{value_type}, vector)

            # vector is a column vector, volume == nothing
            B = Blade{precision_type}(converted_vector)
            @test B isa Blade{precision_type}
            @test B.dim == 3
            @test B.grade == 1
            @test B.basis ≈ reshape(converted_vector / 13,
                                    length(converted_vector), 1)
            @test LinearAlgebra.norm(B.basis) ≈ 1
            @test B.volume ≈ 13

            # vector is a column vector, volume != nothing
            B = Blade{precision_type}(converted_vector, volume=test_volume)
            @test B isa Blade{precision_type}
            @test B.dim == 3
            @test B.grade == 1
            @test B.basis ≈ reshape(converted_vector / 13,
                                    length(converted_vector), 1)
            @test LinearAlgebra.norm(B.basis) ≈ 1
            @test B.volume ≈ test_volume

            # vector is a column vector, volume < default atol
            B = Blade{precision_type}(converted_vector,
                                      volume=blade_atol(precision_type) / 2)
            @test iszero(B)

            # vector is a row vector
            row_vector = reshape(converted_vector, 1, length(converted_vector))
            B = Blade{precision_type}(row_vector)
            @test B.dim == 3
            @test B.grade == 1
            @test B.basis ≈ permutedims(row_vector / 13)
            @test LinearAlgebra.norm(B.basis) ≈ 1
            @test B.volume ≈ 13

            # vector is a zero vector, volume == nothing
            zero_vector = Array{value_type}([0., 0., 0.])
            B = Blade{precision_type}(zero_vector)
            @test iszero(B)

            # vector is a zero vector, volume != nothing
            zero_vector = Array{value_type}([0., 0., 0.])
            B = Blade{precision_type}(zero_vector, volume=test_volume)
            @test iszero(B)

            # norm(vector) < default atol
            small_vector = Vector{value_type}(converted_vector)
            small_vector /= LinearAlgebra.norm(small_vector)
            small_vector *= blade_atol(precision_type) / 2
            B = Blade{precision_type}(small_vector)
            @test iszero(B)

            # norm(vector) < atol
            B = Blade{precision_type}(
                converted_vector,
                atol=LinearAlgebra.norm(converted_vector) + 1)
            @test iszero(B)

            # length(vector) == 0
            B = Blade{precision_type}(rand(0))
            @test iszero(B)
        end

        # --- vector::Vector{<:Integer}
        #
        # Notes
        # -----
        # * We only need to verify the correctness of the type of the result
        #   for a few permutations of the arguments because the computations
        #   for the data fields are independent of the type of the `vectors`
        #   argument. The correctness of the data field computations is
        #   verified in the unit tests for the cases when `vector` has type
        #   Vector{<:AbstractFloat}.

        # subtypes(Signed)
        for value_type in subtypes(Signed)
            # Preparations
            converted_vector = convert(Vector{value_type}, vector)

            # vector is a column vector
            B = Blade{precision_type}(converted_vector)
            @test B isa Blade{precision_type}

            # vector is a column vector, volume < default atol
            B = Blade{precision_type}(converted_vector,
                                      volume=blade_atol(precision_type) / 2)
            @test iszero(B)

            # vector is a row vector
            row_vector = reshape(converted_vector, 1, length(converted_vector))
            B = Blade{precision_type}(row_vector)
            @test B isa Blade{precision_type}

            # vector is a zero vector
            zero_vector = Array{value_type}([0, 0, 0])
            B = Blade{precision_type}(zero_vector)
            @test iszero(B)

            # norm(vector) < default atol
            small_vector = Vector{value_type}(converted_vector)
            small_vector /= LinearAlgebra.norm(small_vector)
            small_vector *= blade_atol(precision_type) / 2
            B = Blade{precision_type}(small_vector)
        end

        # subtypes(Unsigned)
        for value_type in subtypes(Unsigned)
            # Preparations
            converted_vector = convert(Vector{value_type}, abs.(vector))

            # vector is a column vector
            B = Blade{precision_type}(converted_vector)
            @test B isa Blade{precision_type}

            # vector is a column vector, volume < default atol
            B = Blade{precision_type}(converted_vector,
                                      volume=blade_atol(precision_type) / 2)
            @test iszero(B)

            # vector is a row vector
            row_vector = reshape(converted_vector, 1, length(converted_vector))
            B = Blade{precision_type}(row_vector)
            @test B isa Blade{precision_type}

            # vector is a zero vector
            zero_vector = Array{value_type}([0, 0, 0])
            B = Blade{precision_type}(zero_vector)
            @test iszero(B)

            # norm(vector) < default atol
            small_vector = Vector{value_type}(converted_vector)
            small_vector /= LinearAlgebra.norm(small_vector)
            small_vector *= blade_atol(precision_type) / 2
            B = Blade{precision_type}(small_vector)
        end

        # --- Bool

        # Preparations
        converted_vector = Vector{Bool}(vector .!= 0)

        # vector is a column vector
        B = Blade{precision_type}(converted_vector)
        @test B isa Blade{precision_type}

        # vector is a column vector, volume < default atol
        B = Blade{precision_type}(converted_vector,
                                  volume=blade_atol(precision_type) / 2)
        @test iszero(B)

        # vector is a row vector
        row_vector = reshape(converted_vector, 1, length(converted_vector))
        B = Blade{precision_type}(row_vector)
        @test B isa Blade{precision_type}

        # vector is a zero vector
        zero_vector = Array{Bool}([0, 0, 0])
        B = Blade{precision_type}(zero_vector)
        @test iszero(B)
    end
end

@testset "Blade: outer constructors - Blade(vectors::Array)" begin
    #=
      Notes
      -----
      * Test type of constructed instances. Correct construction of instances
        is tested by the inner constructor tests.

      * Test behavior of keyword arguments: `volume`, `atol`, `copy_basis`.
    =#

    # --- Preparations

    vectors = Matrix([3 3; -4 -4; 0 1])
    one_vector = Vector([3; 4; 0])
    pseudoscalar_vectors = Matrix([3 5; 4 12])

    # Select random volume for tests where `volume` != 0
    test_volume = rand() + 1  # add 1 to avoid 0
    test_volume = rand() > 0.5 ? test_volume : -test_volume

    # --- Blade(vectors::Array{T};
    #           volume::Union{Real, Nothing}=nothing,
    #           atol::Real=blade_atol(T)) where {T<:AbstractFloat}

    for precision_type in subtypes(AbstractFloat)
        # --- vectors isa Matrix

        # Preparations
        converted_vectors = convert(Matrix{precision_type}, vectors)

        # norm(blade) > default atol, volume == nothing
        B = Blade(converted_vectors)
        @test B isa Blade{precision_type}

        # norm(blade) > default atol, volume != nothing
        B = Blade(converted_vectors, volume=test_volume)
        @test B isa Blade{precision_type}

        F = LinearAlgebra.qr(converted_vectors)
        signed_norm = prod(LinearAlgebra.diag(F.R))
        @test volume(B) == sign(signed_norm) * precision_type(test_volume)

        # norm(blade) < default atol
        small_blade = blade_atol(precision_type) * converted_vectors / 6
        B = Blade(small_blade)
        @test iszero(B)

        # norm(blade) < atol
        B = Blade(converted_vectors, atol=6)
        @test iszero(B)

        # --- vectors isa Vector

        # Preparations
        converted_one_vector = convert(Vector{precision_type}, one_vector)

        # norm(blade) > default atol, volume == nothing
        B = Blade(converted_one_vector)
        @test B isa Blade{precision_type}

        # norm(blade) > default atol, volume != nothing
        B = Blade(converted_one_vector, volume=test_volume)
        @test B isa Blade{precision_type}
        @test volume(B) == precision_type(test_volume)

        # norm(blade) < default atol
        small_blade = blade_atol(precision_type) * converted_one_vector / 6
        B = Blade(small_blade)
        @test iszero(B)

        # norm(blade) < atol
        B = Blade(converted_one_vector, atol=6)
        @test iszero(B)

        # --- number of vectors == dimension of column space

        # vectors are linearly independent
        converted_vectors = Matrix{precision_type}(pseudoscalar_vectors)
        B = Blade(converted_vectors)
        @test B isa Pseudoscalar{precision_type}

        # vectors are linearly dependent
        converted_vectors = Matrix{precision_type}([2 3 5; 4 6 10; 1 1 2])
        B = Blade(converted_vectors)
        @test iszero(B)
    end

    # --- Blade(vectors::Array{<:Integer};
    #           volume::Union{Real, Nothing}=nothing,
    #           atol::Real=blade_atol(Float64))

    # subtypes(Signed)
    for vectors_type in subtypes(Signed)
        # --- vectors isa Matrix

        # Preparations
        converted_vectors = convert(Matrix{vectors_type}, vectors)

        # norm(blade) > default atol, volume == nothing
        B = Blade(converted_vectors)
        @test B isa Blade{Float64}

        # norm(blade) > default atol, volume != nothing
        B = Blade(converted_vectors, volume=test_volume)
        @test B isa Blade{Float64}

        F = LinearAlgebra.qr(converted_vectors)
        signed_norm = prod(LinearAlgebra.diag(F.R))
        @test volume(B) == sign(signed_norm) * Float64(test_volume)

        # norm(blade) < atol
        B = Blade(converted_vectors, atol=6)
        @test iszero(B)

        # --- vectors isa Vector

        # Preparations
        converted_one_vector = convert(Vector{vectors_type}, one_vector)

        # norm(blade) > default atol, volume == nothing
        B = Blade(converted_one_vector)
        @test B isa Blade{Float64}

        # norm(blade) > default atol, volume != nothing
        B = Blade(converted_one_vector, volume=test_volume)
        @test B isa Blade{Float64}
        @test volume(B) == Float64(test_volume)

        # norm(blade) < atol
        B = Blade(converted_one_vector, atol=6)
        @test iszero(B)

        # --- number of vectors == dimension of column space

        # vectors are linearly independent
        converted_vectors = Matrix{Float64}(pseudoscalar_vectors)
        B = Blade(converted_vectors)
        @test B isa Pseudoscalar{Float64}

        # vectors are linearly dependent
        converted_vectors = Matrix{Float64}([2 3 5; 4 6 10; 1 1 2])
        B = Blade(converted_vectors)
        @test iszero(B)
    end

    # subtypes(Unsigned)
    for vectors_type in subtypes(Unsigned)
        # --- vectors isa Matrix

        # Preparations
        converted_vectors = convert(Matrix{vectors_type}, abs.(vectors))

        # norm(blade) > default atol, volume == nothing
        B = Blade(converted_vectors)
        @test B isa Blade{Float64}

        # norm(blade) > default atol, volume != nothing
        B = Blade(converted_vectors, volume=test_volume)
        @test B isa Blade{Float64}

        F = LinearAlgebra.qr(converted_vectors)
        signed_norm = prod(LinearAlgebra.diag(F.R))
        @test volume(B) == sign(signed_norm) * Float64(test_volume)

        # norm(blade) < atol
        B = Blade(converted_vectors, atol=6)
        @test iszero(B)

        # --- vectors isa Vector

        # Preparations
        converted_one_vector = convert(Vector{vectors_type}, one_vector)

        # norm(blade) > default atol, volume == nothing
        B = Blade(converted_one_vector)
        @test B isa Blade{Float64}

        # norm(blade) > default atol, volume != nothing
        B = Blade(converted_one_vector, volume=test_volume)
        @test B isa Blade{Float64}
        @test volume(B) == Float64(test_volume)

        # norm(blade) < atol
        B = Blade(converted_one_vector, atol=6)
        @test iszero(B)

        # --- number of vectors == dimension of column space

        # vectors are linearly independent
        converted_vectors = Matrix{Float64}(pseudoscalar_vectors)
        B = Blade(converted_vectors)
        @test B isa Pseudoscalar{Float64}

        # vectors are linearly dependent
        converted_vectors = Matrix{Float64}([2 3 5; 4 6 10; 1 1 2])
        B = Blade(converted_vectors)
        @test iszero(B)
    end

    # --- Bool

    # ------ vectors isa Matrix

    # Preparations
    converted_vectors = Matrix{Bool}(vectors .!= 0)

    # norm(blade) > default atol, volume == nothing
    B = Blade(converted_vectors)
    @test B isa Blade{Float64}

    # norm(blade) > default atol, volume != nothing
    B = Blade(converted_vectors, volume=test_volume)
    @test B isa Blade{Float64}

    F = LinearAlgebra.qr(converted_vectors)
    signed_norm = prod(LinearAlgebra.diag(F.R))
    @test volume(B) == sign(signed_norm) * Float64(test_volume)

    # norm(blade) < atol
    B = Blade(converted_vectors, atol=6)
    @test iszero(B)

    # ------ vectors isa Vector

    # Preparations
    converted_one_vector = Vector{Bool}(one_vector .!= 0)

    # norm(blade) > default atol, volume == nothing
    B = Blade(converted_one_vector)
    @test B isa Blade{Float64}

    # norm(blade) > default atol, volume != nothing
    B = Blade(converted_one_vector, volume=test_volume)
    @test B isa Blade{Float64}
    @test volume(B) == Float64(test_volume)

    # norm(blade) < atol
    B = Blade(converted_one_vector, atol=6)
    @test iszero(B)

    # ------ number of vectors == dimension of column space

    # vectors are linearly independent
    converted_vectors = Matrix{Bool}([1 1 0; 0 0 1; 0 1 0])
    B = Blade(converted_vectors)
    @test B isa Pseudoscalar{Float64}
end

@testset "Blade: outer constructors - copy constructors" begin
    #=
      Notes
      -----
      * Test type of constructed instances. Correct construction of instances
        is tested by the inner constructor tests.

      * Test behavior of keyword arguments: `volume`, `atol`, `copy_basis`.
    =#

    # --- Preparations

    vectors = Matrix([3 3; -4 -4; 0 1])
    one_vector = Vector([3; 4; 0])

    # Select random volume for tests where `volume` != 0
    test_volume = rand() + 1  # add 1 to avoid 0
    test_volume = rand() > 0.5 ? test_volume : -test_volume

    # --- Blade{T}(B::Blade{<:AbstractFloat};
    #              volume::Real=volume(B),
    #              atol::Real=blade_atol(T),
    #              copy_basis::Bool=false) where {T<:AbstractFloat}

    for precision_type_src in subtypes(AbstractFloat)
        # Preparations
        vectors = Matrix{precision_type_src}([3 -3; 4 -4; 0 1])
        B = Blade{precision_type_src}(vectors)

        for precision_type_converted in subtypes(AbstractFloat)

            # Construct a Blade representing the same blade as `B`
            B_copy = Blade{precision_type_converted}(B)
            @test B_copy.dim == B.dim
            @test B_copy.grade == B.grade

            if precision_type_src == precision_type_converted
                @test B_copy.basis === B.basis
            else
                @test B_copy.basis ≈ B.basis
                @test B_copy.basis !== B.basis
            end

            @test B_copy.volume ≈ B.volume

            # Construct a Blade representing the same space as `B` with
            # specified volume
            B_copy = Blade{precision_type_converted}(B, volume=test_volume)
            @test B_copy.dim == B.dim
            @test B_copy.grade == B.grade

            if precision_type_src == precision_type_converted
                @test B_copy.basis === B.basis
            else
                @test B_copy.basis ≈ B.basis
                @test B_copy.basis !== B.basis
            end

            @test B_copy.volume == precision_type_converted(test_volume)

            # Construct a Blade representing the same blade as `B` containing
            # a copy of the basis (instead of a reference).
            B_copy = Blade{precision_type_converted}(B, copy_basis=true)
            @test B_copy.dim == B.dim
            @test B_copy.grade == B.grade

            if precision_type_src == precision_type_converted
                @test B_copy.basis == B.basis
            else
                @test B_copy.basis ≈ B.basis
            end
            @test B_copy.basis !== B.basis

            @test B_copy.volume == B.volume
        end
    end

    # --- Blade(B::AbstractBlade{T}; volume=volume(B), copy_basis=false)
    #         where {T<:AbstractFloat}

    for precision_type in subtypes(AbstractFloat)
        # Preparations
        converted_vectors = Matrix{precision_type}(vectors)
        B = Blade(converted_vectors)

        # Construct a Blade representing the same blade as `B`
        B_copy = Blade(B)
        @test B_copy isa Blade{precision_type}
        @test B_copy.basis === B.basis

        # Construct a Blade representing the same blade as `B` with specified
        # volume
        B_copy = Blade(B, volume=test_volume)
        @test B_copy isa Blade{precision_type}
        @test volume(B_copy) == precision_type(test_volume)

        # Construct Blade representing the same blade as `B` containing
        # a copy of the basis (instead of a reference).
        B_copy = Blade(B, copy_basis=true)
        @test B_copy isa Blade{precision_type}
        @test B_copy.basis == B.basis
        @test B_copy.basis !== B.basis
    end
end

@testset "Blade: outer constructors - Scalar constructors" begin
    # --- Preparations

    test_value = rand() + 1  # add 1 to avoid 0
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- Exercise functionality and check results

    @test Blade(test_value) == Scalar(test_value)
    @test Blade(Scalar(test_value)) == Scalar(test_value)
end

# --- Test attribute methods

@testset "Blade: AbstractMultivector attribute functions" begin
    # --- Preparations

    vectors = [3 3; 4 4; 0 1]
    test_dim, test_grade = size(vectors)

    # --- Test basic functions

    for precision_type in subtypes(AbstractFloat)
        # Blade with sign > 0
        B = Blade{precision_type}(vectors)

        @test dim(B) == test_dim
        @test grades(B) == [test_grade]
        @test blades(B) == [B]

        @test norm(B) isa precision_type
        @test norm(B) ≈ 5

        for k in 0:test_dim
            if k == test_grade
                @test B[k] == [B]
            else
                @test B[k] == []
            end
        end
    end
end

@testset "Blade: AbstractBlade attribute functions" begin
    # --- Preparations

    vectors = [3 3; 4 4; 0 1]
    test_grade = size(vectors, 2)

    # --- Test basic functions

    for precision_type in subtypes(AbstractFloat)
        # Blade with sign > 0
        B = Blade{precision_type}(vectors)

        @test grade(B) == test_grade

        F = LinearAlgebra.qr(vectors)
        @test basis(B) ≈ Matrix(F.Q)

        @test volume(B) isa precision_type
        @test volume(B) ≈ 5

        @test sign(B) == 1

        # Blade with sign < 0
        C = Blade(B, volume=-1)
        @test sign(C) == -1
    end
end

# --- Tests for AbstractMultivector interface functions

@testset "Blade: -(B)" begin
    # Preparations
    vectors = Matrix{Float16}([3 3 3 3; 4 4 4 4; 0 1 0 0; 0 0 1 0; 0 0 0 1])
    B = Blade(vectors)

    # B::Blade
    @test -B == Blade(B, volume=-volume(B))
end

@testset "Blade: reverse(B)" begin
    # mod(grade, 4) == 1
    vectors = Vector([3; 4; 0; 0; 0])
    B = Blade(vectors)
    @test reverse(B) === B
    @test B * reverse(B) ≈ norm(B)^2

    # mod(grade, 4) == 2
    vectors = Matrix([3 3; 4 4; 0 1; 0 0; 0 0])
    B = Blade(vectors)
    @test reverse(B) == -B
    @test B * reverse(B) ≈ norm(B)^2

    # mod(grade, 4) == 3
    vectors = Matrix([3 3 3; 4 4 4; 0 1 0; 0 0 1; 0 0 0])
    B = Blade(vectors)
    @test reverse(B) == -B
    @test B * reverse(B) ≈ norm(B)^2

    # mod(grade, 4) == 0
    vectors = Matrix([3 3 3 3; 4 4 4 4; 0 1 0 0; 0 0 1 0; 0 0 0 1])
    B = Blade(vectors)
    @test reverse(B) === B
    @test B * reverse(B) ≈ norm(B)^2
end

@testset "Blade: dual(B)" begin
    # --- Preparations

    # Test values
    test_value = rand()
    test_value = rand() > 0.5 ? test_value : -test_value

    # --- B::Blade

    for dim_B in 10:13
        for grade_B in 2:5
            B = Blade(rand(dim_B, grade_B))

            dual_B = dual(B)

            # --- Check dim, grade, and norm

            @test dim(dual_B) == dim_B
            @test grade(dual_B) == dim_B - grade_B
            @test norm(dual_B) == norm(B)

            # --- Check that B and dual(B) are orthogonal complements

            @test LinearAlgebra.norm(transpose(basis(dual_B)) * basis(B)) <
                10 * eps(Float64)

            # --- Check sign(dual(B))

            # Compute sign of I formed from basis(B) and basis(dual(B))
            sign_Q = sign(LinearAlgebra.det(hcat(basis(B), basis(dual_B))))

            # Compute expected_sign
            expected_sign = sign(B) * sign_Q

            # Account for sign of I^{-1} relative to I
            if mod(dim(B), 4) >= 2
                expected_sign = -expected_sign
            end

            # Account for reversals required to eliminate B
            if mod(grade(B), 4) >= 2
                expected_sign = -expected_sign
            end

            @test sign(dual_B) == expected_sign

            # Check dual(dual_B) = (-1)^(dim_B * (dim_B - 1) / 2) B
            if mod(dim(B), 4) < 2
                @test dual(dual_B) ≈ B
            else
                @test dual(dual_B) ≈ -B
            end
        end
    end
end

@testset "Blade: convert(B)" begin
    # Preparations
    vectors = [3 3; 4 4; 0 1]

    # Tests
    for precision_type_converted in subtypes(AbstractFloat)
        for precision_type_src in subtypes(AbstractFloat)
            # Preparations
            converted_vectors = convert.(precision_type_src, vectors)
            B = Blade{precision_type_src}(converted_vectors)

            # Exercise functionality and check results
            B_converted = convert(precision_type_converted, B)

            @test B_converted isa Blade{precision_type_converted}

            if precision_type_src == precision_type_converted
                @test B_converted === B
            else
                @test B_converted !== B
                @test B_converted ≈ B
            end
        end
    end
end

# --- Tests for AbstractBlade interface functions

@testset "Blade: inv(B)" begin
    # mod(grade, 4) == 1
    vectors = Vector([3; 4; 0; 0; 0])
    B = Blade(vectors)
    expected_inverse = Blade(B, volume=1 / volume(B))
    @test inv(B) ≈ expected_inverse
    @test B * inv(B) ≈ 1

    # mod(grade, 4) == 2
    vectors = Matrix([3 3; 4 4; 0 1; 0 0; 0 0])
    B = Blade(vectors)
    expected_inverse = Blade(B, volume=-1 / volume(B))
    @test inv(B) ≈ expected_inverse
    @test B * inv(B) ≈ 1

    # mod(grade, 4) == 3
    vectors = Matrix([3 3 3; 4 4 4; 0 1 0; 0 0 1; 0 0 0])
    B = Blade(vectors)
    expected_inverse = Blade(B, volume=-1 / volume(B))
    @test inv(B) ≈ expected_inverse
    @test B * inv(B) ≈ 1

    # mod(grade, 4) == 0
    vectors = Matrix([3 3 3 3; 4 4 4 4; 0 1 0 0; 0 0 1 0; 0 0 0 1])
    B = Blade(vectors)
    expected_inverse = Blade(B, volume=1 / volume(B))
    @test inv(B) ≈ expected_inverse
    @test B * inv(B) ≈ 1
end