text:
"""
Creates the covariance matrix by evaluating the kernel function for each pair of passed inputs
symmetric = a parameter stating whether the covariance is guarunteed to be symmetric about the diagonal
"""
function covariance!(
    Σ::DenseArray{T1,2},
    kernel_func::Function,
    x1list::Vector{T1},
    x2list::Vector{T1},
    kernel_hyperparameters::Vector{T1};
    dorder::Vector{T2}=zeros(Int64, 2 + length(kernel_hyperparameters)),
    symmetric::Bool=false
    ) where {T1<:Real, T2<:Integer}

    @assert issorted(x1list)

    # are the list of x's passed identical
    same_x = (x1list == x2list)

    # are the x's passed identical and equally spaced
    if same_x
        spacing = x1list[2:end]-x1list[1:end-1]
        equal_spacing = all((spacing .- spacing[1]) .< (1e-8 * spacing[1]))
    else
        @assert issorted(x2list)
        equal_spacing = false
    end

    x1_length = length(x1list)
    x2_length = length(x2list)

    @assert size(Σ, 1) == x1_length
    @assert size(Σ, 2) == x2_length

    if equal_spacing && symmetric
        # this section is so fast, it isn't worth parallelizing
        kernline = zeros(x1_length)
        for i in 1:x1_length
            kernline[i] = kernel_func(kernel_hyperparameters, x1list[1] - x1list[i], dorder)
        end
        for i in 1:x1_length
            Σ[i, i:end] = kernline[1:(x1_length + 1 - i)]
        end
        return Symmetric(Σ)
    else
        _covariance!(Σ, kernel_func, x1list, x2list, kernel_hyperparameters, dorder, symmetric, same_x, equal_spacing)
    end
end

function _covariance!(
    Σ::SharedArray{T1,2},
    kernel_func::Function,
    x1list::Vector{T1},
    x2list::Vector{T1},
    kernel_hyperparameters::Vector{T1},
    dorder::Vector{T2},
    symmetric::Bool,
    same_x::Bool,
    equal_spacing::Bool
    ) where {T1<:Real, T2<:Integer}

    x1_length = length(x1list)
    x2_length = length(x2list)

    @assert size(Σ) == (x1_length, x2_length)

    if same_x && symmetric
        sendto(workers(), kernel_func=kernel_func, kernel_hyperparameters=kernel_hyperparameters, x1list=x1list, dorder=dorder)
        @sync @distributed for i in 1:length(x1list)
            for j in 1:length(x1list)
                if i <= j; Σ[i, j] = kernel_func(kernel_hyperparameters, x1list[i] - x1list[j], dorder) end
            end
        end
        return Symmetric(Σ)
    else
        sendto(workers(), kernel_func=kernel_func, kernel_hyperparameters=kernel_hyperparameters, x1list=x1list, x2list=x2list, dorder=dorder)
        @sync @distributed for i in 1:length(x1list)
            for j in 1:length(x2list)
                Σ[i, j] = kernel_func(kernel_hyperparameters, x1list[i] - x2list[j], dorder)
            end
        end
        return Σ
    end
end

function _covariance!(
    Σ::Matrix{T1},
    kernel_func::Function,
    x1list::Vector{T1},
    x2list::Vector{T1},
    kernel_hyperparameters::Vector{T1},
    dorder::Vector{T2},
    symmetric::Bool,
    same_x::Bool,
    equal_spacing::Bool
    ) where {T1<:Real, T2<:Integer}

    x1_length = length(x1list)
    x2_length = length(x2list)

    @assert size(Σ) == (x1_length, x2_length)

    if same_x && symmetric
        for i in 1:length(x1list)
            for j in 1:length(x1list)
                if i <= j; Σ[i, j] = kernel_func(kernel_hyperparameters, x1list[i] - x1list[j], dorder) end
            end
        end
        return Symmetric(Σ)
    else
        for i in 1:length(x1list)
            for j in 1:length(x2list)
                Σ[i, j] = kernel_func(kernel_hyperparameters, x1list[i] - x2list[j], dorder)
            end
        end
        return Σ
    end
end

function covariance(
    kernel_func::Function,
    x1list::Vector{T1},
    x2list::Vector{T1},
    kernel_hyperparameters::Vector{T1};
    kwargs...) where {T1<:Real}

    if nworkers()>1
        return covariance!(SharedArray{Float64}(length(x1list), length(x2list)), kernel_func, x1list, x2list, kernel_hyperparameters; kwargs...)
    else
        return covariance!(zeros(length(x1list), length(x2list)), kernel_func, x1list, x2list, kernel_hyperparameters; kwargs...)
    end
    # return covariance!(zeros(length(x1list), length(x2list)), kernel_func, x1list, x2list, kernel_hyperparameters; dorder=dorder, symmetric=symmetric)
end

"""
Calculating the covariance between all outputs for a combination of dependent GPs
written so that the intermediate Σ's don't have to be calculated over and over again
"""
function covariance(
    prob_def::GLO,
    x1list::Vector{<:Real},
    x2list::Vector{<:Real},
    total_hyperparameters::Vector{<:Real};
    dΣdθs_total::Vector{<:Integer}=Int64[],
    kwargs...)

    @assert all(dΣdθs_total .>= 0)
    @assert length(total_hyperparameters) == prob_def.n_kern_hyper + length(prob_def.a0)

    num_coefficients = length(total_hyperparameters) - prob_def.n_kern_hyper
    n_out = prob_def.n_out
    n_dif = prob_def.n_dif
    dΣdθs_kernel = dΣdθs_total .- num_coefficients
    kernel_hyperparameters = total_hyperparameters[(num_coefficients + 1):end]
    # println(length(kernel_hyperparameters))

    x1_length = length(x1list)
    x2_length = length(x2list)

    if nworkers()>1
        holder = SharedArray{Float64}(length(x1list), length(x2list))
    else
        holder = zeros(length(x1list), length(x2list))
    end

    # only calculating each sub-matrix once and using the fact that they should
    # be basically the same if the kernel has been differentiated the same amount of times
    A_list = Vector{AbstractArray{<:Real,2}}(undef, 2 * n_dif - 1)
    dorder = zeros(Int64, 2 + length(kernel_hyperparameters))
    for dΣdθ_kernel in dΣdθs_kernel
        if dΣdθ_kernel > 0; dorder[2 + dΣdθ_kernel] += 1 end
    end
    if !prob_def.kernel_changes_with_output
        for i in 0:(2 * n_dif - 2)
            dorder[1] = rem(i - 1, 2) + 1
            dorder[2] = 2 * div(i - 1, 2)
            # things that have been differentiated an even amount of times are symmetric about t1-t2==0
            if iseven(i)
                A_list[i + 1] = copy(covariance!(holder, prob_def.kernel, x1list, x2list, kernel_hyperparameters; dorder=dorder, symmetric=true))
            else
                A_list[i + 1] = copy(covariance!(holder, prob_def.kernel, x1list, x2list, kernel_hyperparameters; dorder=dorder))
            end
        end
    end
    # assembling the coefficient matrix
    a = reshape(total_hyperparameters[1:num_coefficients], (n_out, n_dif))

    # initializing the multi-output covariance matrix
    Σ = zeros((n_out * x1_length, n_out * x2_length))

    coeff_orders = copy(prob_def.coeff_orders)
    coeff_coeffs = copy(prob_def.coeff_coeffs)
    dif_coefficients!(n_out, n_dif, dΣdθs_total, coeff_orders, coeff_coeffs)
    for i in 1:n_out
        for j in 1:n_out
            # if i <= j  # only need to calculate one set of the off-diagonal blocks
            if prob_def.kernel_changes_with_output
                outputs = [i,j]
                for i in 0:(2 * n_dif - 2)
                    dorder[1] = rem(i - 1, 2) + 1
                    dorder[2] = 2 * div(i - 1, 2)
                    # things that have been differentiated an even amount of times are symmetric about t1-t2==0
                    if iseven(i)
                        A_list[i + 1] = copy(covariance!(holder, (hyper, δ, dord) -> prob_def.kernel(hyper, δ, dord; outputs=outputs), x1list, x2list, kernel_hyperparameters; dorder=dorder, symmetric=true))
                    else
                        A_list[i + 1] = copy(covariance!(holder, (hyper, δ, dord) -> prob_def.kernel(hyper, δ, dord; outputs=outputs), x1list, x2list, kernel_hyperparameters; dorder=dorder))
                    end
                end
            end
                for k in 1:n_dif
                    for l in 1:n_dif
                        # if the coefficient for the GLOM coefficients is non-zero
                        if coeff_coeffs[i, j, k, l] != 0
                            # make it negative or not based on how many times it has been differentiated in the x1 direction
                            A_mat_coeff = coeff_coeffs[i, j, k, l] * powers_of_negative_one(l - 1)
                            for m in 1:n_out
                                for n in 1:n_dif
                                    A_mat_coeff *= a[m, n] ^ coeff_orders[i, j, k, l, m, n]
                                end
                            end
                            # add the properly negative differentiated A matrix from the list
                            Σ[((i - 1) * x1_length + 1):(i * x1_length),
                                ((j - 1) * x2_length + 1):(j * x2_length)] +=
                                A_mat_coeff * A_list[k + l - 1]
                        end
                    end
                end
            # end
        end
    end

    # making the highest variances all along the main diagonals
    #TODO implement in above step (aka at initial assignment)
    Σ_rearranged = zeros((n_out * x1_length, n_out * x2_length))
    for i in 1:x1_length
        for j in 1:x2_length
            Σ_rearranged[1 + (i - 1) * n_out:i * n_out,1 + (j - 1) * n_out:j * n_out] = Σ[i:x1_length:end,:][:, j:x2_length:end]
        end
    end
    # return the symmetrized version of the covariance matrix
    # function corrects for numerical errors and notifies us if our matrix isn't
    # symmetric like it should be
    return symmetric_A(Σ_rearranged; kwargs...)
    # if x1list == x2list
    #     if chol
    #         return ridge_chol(Symmetric(Σ))
    #     else
    #         return Symmetric(Σ)
    #     end
    # else
    #     return Σ
    # end

end

covariance(
    prob_def::GLO,
    total_hyperparameters::Vector{<:Real};
    kwargs...
    ) = covariance(prob_def, prob_def.x_obs, prob_def.x_obs, total_hyperparameters; kwargs...)

"adding measurement noise to Σ_obs"
function Σ_observations(
    Σ::Symmetric{T, Matrix{T}},
    measurement_noise::Vector{T};
    return_both::Bool=false,
    kwargs...) where {T<:Real}

    Σ_obs = symmetric_A(Σ + Diagonal(measurement_noise .* measurement_noise); chol=true, kwargs...)
    return return_both ? (Σ_obs, Σ) : Σ_obs
end

"adding measurement noise to Σ_obs"
function Σ_observations(
    Σ::Symmetric{T, Matrix{T}},
    measurement_covariance::Array{T, 3};
    return_both::Bool=false,
    kwargs...) where {T<:Real}

    n_meas = size(measurement_covariance, 1)
    n_out = size(measurement_covariance, 2)
    @assert n_out == size(measurement_covariance, 3) == size(Σ, 2) / n_meas
    Σ_obs = Matrix(Σ)
    # for i in 1:n_out
    #     for j in 1:n_out
    #         if i <= j; Σ_obs[diagind(Σ_obs, i - 1)[((j - i) * n_meas + 1):((j - i + 1) * n_meas)]] += measurement_covariance[:, i, j] end
    #     end
    # end
    for i in 1:n_meas
        Σ_obs[1 + (i - 1) * n_out : i * n_out, 1 + (i - 1) * n_out : i * n_out] += measurement_covariance[i, :, :]
    end
    Σ_obs = symmetric_A(Σ_obs; chol=true, kwargs...)
    return return_both ? (Σ_obs, Σ) : Σ_obs
end

"adding measurement noise to Σ_obs"
function Σ_observations(
    kernel_func::Function,
    x_obs::Vector{T},
    measurement_noise::Vector{T},
    kernel_hyperparameters::Vector{T};
    ignore_asymmetry::Bool=false,
    return_both::Bool=false
    ) where {T<:Real}

    return Σ_observations(symmetric_A(covariance(kernel_func, x_obs, x_obs, kernel_hyperparameters); ignore_asymmetry=ignore_asymmetry), measurement_noise; return_both=return_both)
end

"adding measurement noise to Σ_obs"
function Σ_observations(
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    ignore_asymmetry::Bool=false,
    return_both::Bool=false
    ) where {T<:Real}

    return prob_def.has_covariance ?
        Σ_observations(symmetric_A(covariance(prob_def, total_hyperparameters); ignore_asymmetry=ignore_asymmetry), prob_def.covariance; return_both=return_both) :
        Σ_observations(symmetric_A(covariance(prob_def, total_hyperparameters); ignore_asymmetry=ignore_asymmetry), prob_def.noise; return_both=return_both)
end

"calculating the standard deviation at each GP posterior point. Algorithm from RW alg. 2.1"
function get_σ(
    L_obs::LowerTriangular{T,Matrix{T}},
    Σ_obs_samp::Union{Transpose{T,Matrix{T}},Symmetric{T,Matrix{T}},Matrix{T}},
    diag_Σ_samp::Vector{T}
    ) where {T<:Real}

    v = L_obs \ Σ_obs_samp
    # thing = diag_Σ_samp - [dot(v[:, i], v[:, i]) for i in 1:length(diag_Σ_samp)]
    # thing[abs.(thing) .< 1e-10] .= 0
    # return sqrt.(thing)  # σ
    return sqrt.(diag_Σ_samp - [dot(v[:, i], v[:, i]) for i in 1:length(diag_Σ_samp)])  # σ
    # return sqrt.(diag_Σ_samp - diag(Σ_samp_obs * (Σ_obs \ Σ_obs_samp)))  # much slower
end

function get_σ(
    prob_def::GLO,
    x_samp::Vector{T},
    total_hyperparameters::Vector{T}
    ) where {T<:Real}

    (Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp) = covariance_permutations(prob_def, x_samp, total_hyperparameters)
    return get_σ(ridge_chol(Σ_obs).L, Σ_obs_samp, diag(Σ_samp))
end

"calcuate all of the different versions of the covariance matrices for measured and sampled points"
function covariance_permutations(
    kernel_func::Function,
    x_obs::Vector{T},
    x_samp::Vector{T},
    measurement_noise::Vector{T},
    kernel_hyperparameters::Vector{T};
    return_both::Bool=false
    ) where {T<:Real}

    Σ_samp = covariance(kernel_func, x_samp, x_samp, kernel_hyperparameters)
    Σ_samp_obs = covariance(kernel_func, x_samp, x_obs, kernel_hyperparameters; ignore_asymmetry=true)
    Σ_obs_samp = transpose(Σ_samp_obs)
    if return_both
        Σ_obs, Σ_obs_raw = Σ_observations(kernel_func, x_obs, measurement_noise, kernel_hyperparameters; return_both=return_both)
        return Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp, Σ_obs_raw
    else
        Σ_obs = Σ_observations(kernel_func, x_obs, measurement_noise, kernel_hyperparameters; return_both=return_both)
        return Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp
    end
end

"calcuate all of the different versions of the covariance matrices for measured and sampled points"
function covariance_permutations(
    prob_def::GLO,
    x_samp::Vector{T},
    total_hyperparameters::Vector{T};
    return_both::Bool=false
    ) where {T<:Real}

    Σ_samp = covariance(prob_def, x_samp, x_samp, total_hyperparameters)
    Σ_samp_obs = covariance(prob_def, x_samp, prob_def.x_obs, total_hyperparameters; ignore_asymmetry=true)
    Σ_obs_samp = transpose(Σ_samp_obs)
    if return_both
        Σ_obs, Σ_obs_raw = Σ_observations(prob_def, total_hyperparameters; return_both=return_both)
        return Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp, Σ_obs_raw
    else
        Σ_obs = Σ_observations(prob_def, total_hyperparameters; return_both=return_both)
        return Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp
    end
end

"Condition the GP on data"
function GP_posteriors_from_covariances(
    y_obs::Vector{T},
    Σ_samp::Union{Cholesky{T,Matrix{T}},Symmetric{T,Matrix{T}},Matrix{T}},
    Σ_obs::Cholesky{T,Matrix{T}},
    Σ_samp_obs::Union{Symmetric{T,Matrix{T}},Matrix{T}},
    Σ_obs_samp::Union{Transpose{T,Matrix{T}},Symmetric{T,Matrix{T}},Matrix{T}},
    Σ_obs_raw::Symmetric{T,Matrix{T}};
    return_Σ::Bool=true,
    kwargs...
    ) where {T<:Real}

    # posterior mean calcuation from RW alg. 2.1

    # these are all equivalent but have different computational costs
    # α = inv(Σ_obs) * y_obs
    # α = transpose(L) \ (L \ y_obs)
    α = Σ_obs \ y_obs

    mean_post = Σ_samp_obs * α
    mean_post_obs = Σ_obs_raw * α

    # posterior standard deviation calcuation from RW alg. 2.1
    σ = get_σ(Σ_obs.L, Σ_obs_samp, diag(Σ_samp))

    # posterior covariance calculation is from eq. 2.24 of RW
    if return_Σ
        Σ_post = symmetric_A(Σ_samp - (Σ_samp_obs * (Σ_obs \ Σ_obs_samp)); kwargs...)
        return mean_post, σ, mean_post_obs, Σ_post
    else
        return mean_post, σ, mean_post_obs
    end

end

"Condition the GP on data"
function GP_posteriors_from_covariances(
    y_obs::Vector{T},
    Σ_samp::Union{Cholesky{T,Matrix{T}},Symmetric{T,Matrix{T}},Matrix{T}},
    Σ_obs::Cholesky{T,Matrix{T}},
    Σ_samp_obs::Union{Symmetric{T,Matrix{T}},Matrix{T}},
    Σ_obs_samp::Union{Transpose{T,Matrix{T}},Symmetric{T,Matrix{T}},Matrix{T}};
    return_Σ::Bool=true,
    kwargs...
    ) where {T<:Real}

    # posterior mean calcuation from RW alg. 2.1

    # these are all equivalent but have different computational costs
    # α = inv(Σ_obs) * y_obs
    # α = transpose(L) \ (L \ y_obs)
    α = Σ_obs \ y_obs

    mean_post = Σ_samp_obs * α

    # posterior standard deviation calcuation from RW alg. 2.1
    σ = get_σ(Σ_obs.L, Σ_obs_samp, diag(Σ_samp))

    # posterior covariance calculation is from eq. 2.24 of RW
    if return_Σ
        Σ_post = symmetric_A(Σ_samp - (Σ_samp_obs * (Σ_obs \ Σ_obs_samp)); kwargs...)
        return mean_post, σ, Σ_post
    else
        return mean_post, σ
    end

end

function GP_posteriors(
    kernel_func::Function,
    x_obs::Vector{T},
    y_obs::Vector{T},
    x_samp::Vector{T},
    measurement_noise::Vector{T},
    total_hyperparameters::Vector{T};
    return_mean_obs::Bool=false,
    kwargs...) where {T<:Real}

    if return_mean_obs
        (Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp, Σ_obs_raw) = covariance_permutations(kernel_func, x_obs, x_samp, measurement_noise, total_hyperparameters; return_both=return_mean_obs)
        return GP_posteriors_from_covariances(y_obs, Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp, Σ_obs_raw; kwargs...)
    else
        (Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp) = covariance_permutations(kernel_func, x_obs, x_samp, measurement_noise, total_hyperparameters; return_both=return_mean_obs)
        return GP_posteriors_from_covariances(y_obs, Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp; kwargs...)
    end
end

function GP_posteriors(
    prob_def::GLO,
    x_samp::Vector{T},
    total_hyperparameters::Vector{T};
    return_mean_obs::Bool=false,
    y_obs::Vector{T}=prob_def.y_obs,
    kwargs...) where {T<:Real}

    if return_mean_obs
        (Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp, Σ_obs_raw) = covariance_permutations(prob_def, x_samp, total_hyperparameters; return_both=return_mean_obs)
        return GP_posteriors_from_covariances(y_obs, Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp, Σ_obs_raw; kwargs...)
    else
        (Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp) = covariance_permutations(prob_def, x_samp, total_hyperparameters; return_both=return_mean_obs)
        return GP_posteriors_from_covariances(y_obs, Σ_samp, Σ_obs, Σ_samp_obs, Σ_obs_samp; kwargs...)
    end

end

"""
Find the powers that each GLOM coefficient is taken to for each part of the
matrix construction before differentiating by any hyperparameters.

Parameters:

n_out (int): Amount of dimensions being fit
n_dif (int): Amount of GP time derivatives are in the GLOM model being used
a (matrix): The coefficients for the GLOM model

Returns:
6D matrix: Filled with integers for what power each coefficient is taken to in
    the construction of a given block of the total covariance matrix.
    For example, coeff_orders[1,1,2,3,:,:] would tell you the powers of each
    coefficient (:,:) that are multiplied by the covariance matrix constructed
    by evaluating the partial derivative of the kernel (once by t1 and twice by
    t2) at every pair of time points (2,3) in the construction of the first
    block of the total covariance matrix (1,1)
4D matrix: Filled with ones anywhere that coeff_orders indicates that
    coefficients exists to multiply a given covariance matrix for a given block

"""
function coefficient_orders(
    n_out::Integer,
    n_dif::Integer;
    a::Matrix{T}=ones(n_out, n_dif)
    ) where {T<:Real}

    @assert size(a) == (n_out, n_dif)

    # ((output pair), (which A matrix to use), (which a coefficent to use))
    coeff_orders = zeros(Int64, n_out, n_out, n_dif, n_dif, n_out, n_dif)
    coeff_coeffs = zeros(Int64, n_out, n_out, n_dif, n_dif)
    for i in 1:n_out
        for j in 1:n_out
            for k in 1:n_dif
                for l in 1:n_dif
                    for m in 1:n_out
                        for n in 1:n_dif
                            if a[m, n] != 0
                                if [m, n] == [i, k]; coeff_orders[i, j, k, l, m, n] += 1 end
                                if [m, n] == [j, l]; coeff_orders[i, j, k, l, m, n] += 1 end
                            end
                        end
                    end
                    # There should be two coefficients being applied to every
                    # matrix. If there are less than two, that means one of the
                    # coefficients was zero, so we should just set them both to
                    # zero
                    if sum(coeff_orders[i, j, k, l, :, :]) != 2
                        coeff_orders[i, j, k, l, :, :] .= 0
                    else
                        coeff_coeffs[i, j, k, l] = 1
                    end
                end
            end
        end
    end

    return coeff_orders, coeff_coeffs

end

"""
Getting the coefficients for constructing differentiated versions of the kernel
using the powers that each coefficient is taken to for each part of the matrix
construction

Parameters:

n_out (int): Amount of dimensions being fit
n_dif (int): Amount of GP time derivatives are in the GLOM model being used
dΣdθ_total (matrix): The coefficients for the GLOM model
coeff_orders (6D matrix): Filled with integers for what power each coefficient
    is taken to in the construction of a given block of the total covariance
    matrix. See coefficient_orders()
coeff_coeffs (4D matrix): Filled with ones and zeros based on which time
    differentiated covariance matrices are added to which blocks

Returns:
Only modifies the passed coeff_orders and coeff_coeffs matrices

"""
function dif_coefficients!(
    n_out::Integer,
    n_dif::Integer,
    dΣdθ_total::Integer,
    coeff_orders::AbstractArray{T,6},
    coeff_coeffs::AbstractArray{T,4}
    ) where {T<:Integer}

    @assert size(coeff_orders) == (n_out, n_out, n_dif, n_dif, n_out, n_dif)
    @assert size(coeff_orders) == (n_out, n_out, n_dif, n_dif, n_out, n_dif)

    # only do something if a derivative is being taken
    if dΣdθ_total > 0 && dΣdθ_total <= (n_out*n_dif)
        proper_indices = [((dΣdθ_total - 1) % n_out) + 1, div(dΣdθ_total - 1, n_out) + 1]
        for i in 1:n_out
            for j in 1:n_out
                for k in 1:n_dif
                    for l in 1:n_dif
                        if coeff_orders[i, j, k, l, proper_indices[1], proper_indices[2]] != 0
                            coeff_coeffs[i, j, k, l] *= coeff_orders[i, j, k, l, proper_indices[1], proper_indices[2]]
                            coeff_orders[i, j, k, l, proper_indices[1], proper_indices[2]] -= 1
                        else
                            coeff_coeffs[i, j, k, l] = 0
                            coeff_orders[i, j, k, l, :, :] .= 0
                        end
                    end
                end
            end
        end
    end
end

function dif_coefficients!(
    n_out::Integer,
    n_dif::Integer,
    dΣdθs_total::Vector{T},
    coeff_orders::AbstractArray{T,6},
    coeff_coeffs::AbstractArray{T,4}
    ) where {T<:Integer}

    for dΣdθ_total in dΣdθs_total
        dif_coefficients!(n_out, n_dif, dΣdθ_total, coeff_orders, coeff_coeffs)
    end
end

# TODO: rename nlogL_normalization
"""
Multivariate normal negative log  likelihood

Parameters:
Σ (Cholesky factorized object): The covariance matrix
y (vector): The values
α (vector): inv(Σ_obs) * y_obs

Returns:
float: the negative log marginal likelihood
"""
function nlogL(
    Σ::Cholesky{T,Matrix{T}},
    y::Vector{T};
    α::Vector{T} = Σ \ y,
    nlogL_normalization::T=logdet(Σ)+length(y)*log(2*π)
    ) where {T<:Real}

    # n = length(y)

    # 2 times negative goodness of fit term
    data_fit = transpose(y) * α
    # 2 times negative complexity penalization term
    # complexity_penalty = log(det(Σ_obs))
    # complexity_penalty = logdet(Σ_obs)  # half memory but twice the time
    # 2 times negative normalization term (doesn't affect fitting)
    # normalization = n * log(2 * π)

    return (data_fit + nlogL_normalization) / 2

end

"""
First partial derivative of the GP negative log marginal likelihood w.r.t. GP
hyperparameters
(see eq. 5.9 in Rasmussen and Williams 2006)

Parameters:
y (vector): The observations at each time point
α (vector): inv(Σ_obs) * y
β (matrix): inv(Σ_obs) * dΣ_dθ where dΣ_dθ is the partial derivative of the
    covariance matrix Σ_obs w.r.t. a hyperparameter

Returns:
float: the partial derivative of the negative log marginal likelihood w.r.t. the
    hyperparameter used in the calculation of β

"""
function dnlogLdθ(
    y::Vector{T},
    α::Vector{T},
    β::Matrix{T}
    ) where {T<:Real}

    # 2 times negative derivative of goodness of fit term
    data_fit = -(transpose(y) * β * α)
    # 2 times negative derivative of complexity penalization term
    complexity_penalty = tr(β)

    # return -1 / 2 * tr((α * transpose(α) - inv(Σ_obs)) * dΣ_dθj)
    return (data_fit + complexity_penalty) / 2

end

"""
First partial derivative of the GP negative log marginal likelihood w.r.t.
parameters that affect y

Parameters:
y1 (vector): The derivative of observations at each time point
α (vector): inv(Σ_obs) * y

Returns:
float: the partial derivative of the negative log marginal likelihood w.r.t. the
    parameter used in the calculation of y1

"""
function dnlogLdθ(
    y1::Vector{T},
    α::Vector{T}
    ) where {T<:Real}

    return transpose(y1) * α
    # return transpose(y) * α1
    # return (transpose(y1) * Σ^-1 * y + transpose(y) * Σ^-1 * y1) / 2
end

"Returns gradient of nlogL in GP hyperparameters"
∇nlogL(y::Vector{T}, α::Vector{T}, βs::Vector{Matrix{T}}) where {T<:Real} = [dnlogLdθ(y, α, β) for β in βs]

"""
Second partial derivative of the GP negative log marginal likelihood w.r.t. two
GP hyperparameters.
Calculated with help from rules found on page 7 of the matrix cookbook
(https://www.ics.uci.edu/~welling/teaching/KernelsICS273B/MatrixCookBook.pdf)

Parameters:
y (vector): The observations at each time point
α (vector): inv(Σ_obs) * y
β1 (matrix): inv(Σ_obs) * dΣ_dθ1 where dΣ_dθ1 is the partial derivative of the
    covariance matrix Σ_obs w.r.t. a hyperparameter
β2 (matrix): inv(Σ_obs) * dΣ_dθ2 where dΣ_dθ2 is the partial derivative of the
    covariance matrix Σ_obs w.r.t. another hyperparameter
β12 (matrix): inv(Σ_obs) * d2Σ_dθ1dθ2 where d2Σ_dθ1dθ2 is the partial
    derivative of the covariance matrix Σ_obs w.r.t. both of the hyperparameters
    being considered

Returns:
float: the partial derivative of the negative log marginal likelihood w.r.t. the
    hyperparameters used in the calculation of the inputs

"""
function d2nlogLdθ(
    y::Vector{T},
    α::Vector{T},
    β1::Matrix{T},
    β2::Matrix{T},
    β12::Matrix{T}
    ) where {T<:Real}

    β12mβ2β1 = β12 - β2 * β1

    # 2 times negative second derivative of goodness of fit term
    data_fit = -(transpose(y) * (β12mβ2β1 - β1 * β2) * α)
    # 2 times negative derivative of complexity penalization term
    complexity_penalty = tr(β12mβ2β1)

    return (data_fit + complexity_penalty) / 2

end

"""
Second partial derivative of the GP negative log marginal likelihood w.r.t. two
parameters that affect y.
Calculated with help from rules found on page 7 of the matrix cookbook
(https://www.ics.uci.edu/~welling/teaching/KernelsICS273B/MatrixCookBook.pdf)

Parameters:
y2 (vector): The derivative of observations at each time point w.r.t. the second
    relevant parameter
y12 (vector): The derivative of observations at each time point w.r.t. the both
    relevant parameters
α (vector): inv(Σ_obs) * y
α1 (vector): inv(Σ_obs) * y1

Returns:
float: the partial derivative of the negative log marginal likelihood w.r.t. the
    hyperparameters used in the calculation of the inputs

"""
function d2nlogLdθ(
    y2::Vector{T},
    y12::Vector{T},
    α::Vector{T},
    α1::Vector{T}
    ) where {T<:Real}

    return transpose(y12) * α + transpose(y2) * α1
    # return (transpose(y12) * α + transpose(y1) * α2 + transpose(y2) * α1 + transpose(y) * α12) / 2

end

"""
Second partial derivative of the GP negative log marginal likelihood w.r.t. a
GP hyperparameter and a parameter thats affect y.
Calculated with help from rules found on page 7 of the matrix cookbook
(https://www.ics.uci.edu/~welling/teaching/KernelsICS273B/MatrixCookBook.pdf)

Parameters:
y (vector): The observations at each time point
y1 (vector): The derivative of observations at each time point w.r.t the
    parameter that affects y
α (vector): inv(Σ_obs) * y
α1 (vector): inv(Σ_obs) * y1
β2 (matrix): inv(Σ_obs) * dΣ_dθ2 where dΣ_dθ2 is the partial derivative of the
    covariance matrix Σ_obs w.r.t. the desired GP hyperparameter

Returns:
float: the partial derivative of the negative log marginal likelihood w.r.t. the
    hyperparameters used in the calculation of the inputs

"""
function d2nlogLdθ(
    y::Vector{T},
    y1::Vector{T},
    α::Vector{T},
    α1::Vector{T},
    β2::Matrix{T}
    ) where {T<:Real}

    return -(transpose(y1) * β2 * α + transpose(y) * β2 * α1) / 2

end

# struct nlogL_matrix_workspace
#     nlogL_hyperparameters::Vector{<:Real}
#     Σ_obs::Cholesky{T,Matrix{T}} where T<:Real
#     y_obs::Vector{<:Real}
#     α::Vector{<:Real}
#     ∇nlogL_hyperparameters::Vector{<:Real}
#     βs::Vector{Matrix{<:Real}}
struct nlogL_matrix_workspace{T<:Real}
    nlogL_hyperparameters::Vector{T}
    Σ_obs::Cholesky{T,Matrix{T}}
    ∇nlogL_hyperparameters::Vector{T}
    βs::Vector{Matrix{T}}

    function nlogL_matrix_workspace(prob_def::GLO,
        total_hyperparameters::Vector{<:Real})

        total_hyper = reconstruct_total_hyperparameters(prob_def::GLO, total_hyperparameters)
        Σ_obs = Σ_observations(prob_def, total_hyper; ignore_asymmetry=true)
        return nlogL_matrix_workspace(
            total_hyper,
            Σ_obs,
            copy(total_hyper),
            [Σ_obs \ covariance(prob_def, total_hyper; dΣdθs_total=[i]) for i in prob_def.non_zero_hyper_inds])
    end
    nlogL_matrix_workspace(
        nlogL_hyperparameters::Vector{T},
        Σ_obs::Cholesky{T,Matrix{T}},
        ∇nlogL_hyperparameters::Vector{T},
        βs::Vector{Matrix{T}}
        ) where T<:Real = new{typeof(nlogL_hyperparameters[1])}(nlogL_hyperparameters, Σ_obs, ∇nlogL_hyperparameters, βs)
end

# "Ensure that the passed matrix workspace parameters are what we expect them to be"
# function check_matrix_workspace(
#     nlogL_hyperparameters::Vector{T}
#     Σ_obs::Cholesky{T,Matrix{T}}
#     y_obs::Vector{T}
#     α::Vector{T}
#     ∇nlogL_hyperparameters::Vector{T}
#     βs::Vector{Matrix{T}}
#     ) where {T<:Real}
#
#     @assert length(y_obs) == length(α) == size(Σ_obs, 1) == size(Σ_obs, 2) = size(βs[1], 1) == size(βs[1], 2)
#     @assert length(nlogL_hyperparameters) == length(∇nlogL_hyperparameters)
#     @assert length(findall(!iszero, nlogL_hyperparameters)) == length(findall(!iszero, ∇nlogL_hyperparameters)) == size(βs, 1)
# end

# function init_nlogL_matrix_workspace(
#     prob_def::GLO,
#     total_hyperparameters::Vector{<:Real})
#
#     total_hyper = reconstruct_total_hyperparameters(prob_def, total_hyperparameters)
#     Σ_obs = Σ_observations(prob_def, total_hyper; ignore_asymmetry=true)
#     return nlogL_matrix_workspace(
#         total_hyper,
#         Σ_obs,
#         copy(prob_def.y_obs),
#         Σ_obs \ prob_def.y_obs,
#         copy(total_hyper),
#         [Σ_obs \ covariance(prob_def, total_hyper; dΣdθs_total=[i]) for i in findall(!iszero, total_hyper)])
# end

"Calculates the quantities shared by the nlogL and ∇nlogL calculations"
function calculate_shared_nlogL_matrices(
    prob_def::GLO,
    non_zero_hyperparameters::Vector{<:Real};
    Σ_obs::Cholesky{T,Matrix{T}} where T<:Real=Σ_observations(prob_def, reconstruct_total_hyperparameters(prob_def, non_zero_hyperparameters); ignore_asymmetry=true))

    # this allows us to prevent the optimizer from seeing the constant zero coefficients
    total_hyperparameters = reconstruct_total_hyperparameters(prob_def, non_zero_hyperparameters)

    return total_hyperparameters, Σ_obs

end

function calculate_shared_nlogL_matrices!(
    workspace::nlogL_matrix_workspace,
    prob_def::GLO,
    non_zero_hyperparameters::Vector{<:Real})

    # this allows us to prevent the optimizer from seeing the constant zero coefficients
    total_hyperparameters = reconstruct_total_hyperparameters(prob_def, non_zero_hyperparameters)

    if workspace.nlogL_hyperparameters != total_hyperparameters
        workspace.nlogL_hyperparameters[:] = total_hyperparameters
        workspace.Σ_obs.factors[:,:] = Σ_observations(prob_def, total_hyperparameters).factors
    end

end

"Calculates the quantities shared by the ∇nlogL and ∇∇nlogL calculations"
function calculate_shared_∇nlogL_matrices(
    prob_def::GLO,
    non_zero_hyperparameters::Vector{<:Real};
    kwargs...)

    total_hyperparameters, Σ_obs = calculate_shared_nlogL_matrices(prob_def, non_zero_hyperparameters; kwargs...)

    βs = [Σ_obs \ covariance(prob_def, total_hyperparameters; dΣdθs_total=[i]) for i in prob_def.non_zero_hyper_inds]

    return total_hyperparameters, Σ_obs, βs

end

function calculate_shared_∇nlogL_matrices!(
    workspace::nlogL_matrix_workspace,
    prob_def::GLO,
    non_zero_hyperparameters::Vector{<:Real})

    calculate_shared_nlogL_matrices!(workspace, prob_def, non_zero_hyperparameters)

    if workspace.∇nlogL_hyperparameters != workspace.nlogL_hyperparameters
        workspace.∇nlogL_hyperparameters[:] = workspace.nlogL_hyperparameters
        workspace.βs[:] = [workspace.Σ_obs \ covariance(prob_def, workspace.∇nlogL_hyperparameters; dΣdθs_total=[i]) for i in findall(!iszero, workspace.∇nlogL_hyperparameters)]
    end

end

"""
Tries to include the specified kernel, assuming it was included with GLOM
Returns the kernel function and number of hyperparameters it uses
"""
function include_kernel(kernel_name::AbstractString)
    if !occursin("_kernel", kernel_name)
        kernel_name *= "_kernel"
    end
    return include(joinpath(pathof(GPLinearODEMaker),"..","kernels", kernel_name * ".jl"))
end

"""
Make it easy to run the covariance calculations on many processors
Makes sure every worker has access to kernel function
"""
function prep_parallel_covariance(
    kernel_name::AbstractString,
    kernel_path::AbstractString;
    kwargs...)

    prep_parallel(; kwargs...)
    sendto(workers(), kernel_name=kernel_name)
    @everywhere include(kernel_path)
end

"""
Make it easy to run the covariance calculations on many processors
Makes sure every worker has access to kernel function
"""
function prep_parallel_covariance(
    kernel_name::AbstractString;
    kwargs...)

    prep_parallel(; kwargs...)
    sendto(workers(), kernel_name=kernel_name)
    @everywhere include_kernel(kernel_name)
end

"nlogL for GLOM GP"
function nlogL_GLOM(
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    y_obs::Vector{T}=copy(prob_def.y_obs),
    kwargs...) where {T<:Real}

    total_hyperparameters, Σ_obs = calculate_shared_nlogL_matrices(
        prob_def, remove_zeros(total_hyperparameters); kwargs...)

    return nlogL(Σ_obs, y_obs)
end

"nlogL for GLOM GP"
function nlogL_GLOM!(
    workspace::nlogL_matrix_workspace,
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    y_obs::Vector{T}=copy(prob_def.y_obs)
    ) where {T<:Real}

    calculate_shared_nlogL_matrices!(workspace, prob_def, total_hyperparameters)
    return nlogL(workspace.Σ_obs, y_obs)
end

"Returns gradient of nlogL for non-zero hyperparameters"
function ∇nlogL_GLOM(
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    y_obs::Vector{T}=copy(prob_def.y_obs),
    kwargs...) where {T<:Real}

    total_hyperparameters, Σ_obs, βs = calculate_shared_∇nlogL_matrices(
        prob_def, remove_zeros(total_hyperparameters); kwargs...)

    return ∇nlogL(y_obs, Σ_obs \ y_obs, βs)

end

"Returns gradient of nlogL for non-zero hyperparameters"
function ∇nlogL_GLOM!(
    workspace::nlogL_matrix_workspace,
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    y_obs::Vector{T}=copy(prob_def.y_obs)
    ) where {T<:Real}

    calculate_shared_∇nlogL_matrices!(workspace, prob_def, total_hyperparameters)

    return ∇nlogL(y_obs, workspace.Σ_obs \ y_obs, workspace.βs)

end

"Replaces H with Hessian of nlogL for non-zero hyperparameters"
function ∇∇nlogL_GLOM(
    prob_def::GLO,
    total_hyperparameters::Vector{T},
    Σ_obs::Cholesky{T,Matrix{T}},
    y_obs::Vector{T},
    α::Vector{T},
    βs::Array{Matrix{T},1}
    ) where {T<:Real}

    non_zero_inds = copy(prob_def.non_zero_hyper_inds)
    H = zeros(length(non_zero_inds), length(non_zero_inds))

    for (i, nzind1) in enumerate(non_zero_inds)
        for (j, nzind2) in enumerate(non_zero_inds)
            if i <= j
                H[i, j] = d2nlogLdθ(y_obs, α, βs[i], βs[j], Σ_obs \ covariance(prob_def, total_hyperparameters; dΣdθs_total=[nzind1, nzind2]))
            end
        end
    end

    H = Symmetric(H)

end

"Replaces H with Hessian of nlogL for non-zero hyperparameters"
function ∇∇nlogL_GLOM(
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    y_obs::Vector{T}=copy(prob_def.y_obs),
    kwargs...) where {T<:Real}

    total_hyperparameters, Σ_obs, βs = calculate_shared_∇nlogL_matrices(
        prob_def, remove_zeros(total_hyperparameters); kwargs...)

    return ∇∇nlogL_GLOM(prob_def, total_hyperparameters, Σ_obs, y_obs, Σ_obs \ y_obs, βs)

end

"Returns gradient of nlogL for non-zero hyperparameters"
function ∇∇nlogL_GLOM!(
    workspace::nlogL_matrix_workspace,
    prob_def::GLO,
    total_hyperparameters::Vector{T};
    y_obs::Vector{T}=copy(prob_def.y_obs)
    ) where {T<:Real}

    calculate_shared_∇nlogL_matrices!(workspace, prob_def, total_hyperparameters)

    return ∇∇nlogL_GLOM(prob_def, workspace.nlogL_hyperparameters, workspace.Σ_obs, y_obs, workspace.Σ_obs \ y_obs, workspace.βs)

end

"reinsert the zero coefficients into the non-zero hyperparameter list if needed"
function reconstruct_total_hyperparameters(
    prob_def::GLO,
    hyperparameters::Vector{T}
    ) where {T<:Real}

    if length(hyperparameters)!=(prob_def.n_kern_hyper + length(prob_def.a0))
        new_coeff_array = reconstruct_array(hyperparameters[1:end - prob_def.n_kern_hyper], prob_def.a0)
        coefficient_hyperparameters = collect(Iterators.flatten(new_coeff_array))
        total_hyperparameters = append!(coefficient_hyperparameters, hyperparameters[end - prob_def.n_kern_hyper + 1:end])
    else
        total_hyperparameters = copy(hyperparameters)
    end

    @assert length(total_hyperparameters)==(prob_def.n_kern_hyper + length(prob_def.a0))

    return total_hyperparameters

end


text:
function bollingerbands{T,N}(ta::TimeArray{T,N}, ma::Int, width::Float64)
    tama   = sma(ta, ma)
    upband = tama .+ moving(ta, std, ma) .* width .* sqrt((ma-1)/ma) # take out Bessel correction, per algorithm
    dnband = tama .- moving(ta, std, ma) .* width .* sqrt((ma-1)/ma)
    bands  =  merge(upband, dnband)
    merge(bands, tama, col_names = ["up", "down", "mean"])
end

bollingerbands{T,N}(ta::TimeArray{T,N}) = bollingerbands(ta, 20, 2.0)

function truerange{T,N}(ohlc::TimeArray{T,N}; h="High", l="Low", c="Close")
    highs    = merge(ohlc[h], lag(ohlc[c]))
    lows     = merge(ohlc[l], lag(ohlc[c]))
    truehigh = TimeArray(highs.timestamp, maximum(highs.values, 2), ["hi"], highs.meta)
    truelow  = TimeArray(lows.timestamp,  minimum(lows.values, 2),  ["lo"], lows.meta)
    truehigh .- truelow
end

function atr{T,N}(ohlc::TimeArray{T,N}, n::Int; h="High", l="Low", c="Close")
    # atr was inveted by Wilder, so only his ema is currently supported
    res = ema(truerange(ohlc), n, wilder=true)
    TimeArray(res.timestamp, res.values, ["atr"], ohlc.meta)
end

atr{T,N}(ta::TimeArray{T,N}) = atr(ta, 14)

# function keltnerbands{T,N}(ohlc::TimeArray{T,N}, n::Int)
# 	typ = typical(ohlc)
# 	rng = ohlc["High"] .- ohlc["Low"]
# 	rma = sma(rng, n)
#
# 	kma     = sma(typ, n)
# 	tstamps = kma.timestamp
#
# 	kma = TimeArray(tstamps, kma.values, ["kma"])
# 	kup = TimeArray(tstamps, (kma.+rma).values, ["kup"])
# 	kdn = TimeArray(tstamps, (kma.-rma).values, ["kdn"])
#
# 	merge(kma, merge(kup, kdn))
# end

keltnerbands{T,N}(ohlc::TimeArray{T,N}) = keltnerbands(ohlc, 10)

# # function chaikinvolatility{T,N}(ta::TimeArray{T,N}, n::Int)
# #   #code here
# # end

ext: jl


text:
# Autogenerated wrapper script for prrte_jll for x86_64-linux-musl
export prte, prte_info, prterun, prun

using libevent_jll
using Hwloc_jll
using PMIx_jll
JLLWrappers.@generate_wrapper_header("prrte")
JLLWrappers.@declare_executable_product(prte)
JLLWrappers.@declare_executable_product(prte_info)
JLLWrappers.@declare_executable_product(prterun)
JLLWrappers.@declare_executable_product(prun)
function __init__()
    JLLWrappers.@generate_init_header(libevent_jll, Hwloc_jll, PMIx_jll)
    JLLWrappers.@init_executable_product(
        prte,
        "bin/prte",
    )

    JLLWrappers.@init_executable_product(
        prte_info,
        "bin/prte_info",
    )

    JLLWrappers.@init_executable_product(
        prterun,
        "bin/prterun",
    )

    JLLWrappers.@init_executable_product(
        prun,
        "bin/prun",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using PiGPIOC

function send_16(spi, data::AbstractArray{UInt8})
    cstr = Cstring(pointer(data))
    ret = PiGPIOC.spiWrite(spi, cstr, 2)
    if ret < 0
        throw("spiWrite returns error with status: $ret")
    end

    nothing
end

function clear_all(spi)
    send_16(spi, [0b0000_0001, 0b0000_0000])
    send_16(spi, [0b0000_0010, 0b0000_0000])
    send_16(spi, [0b0000_0011, 0b0000_0000])
    send_16(spi, [0b0000_0100, 0b0000_0000])
    send_16(spi, [0b0000_0101, 0b0000_0000])
    send_16(spi, [0b0000_0110, 0b0000_0000])
    send_16(spi, [0b0000_0111, 0b0000_0000])
    send_16(spi, [0b0000_1000, 0b0000_0000])
end

println("Starting...")

if PiGPIOC.gpioInitialise() < 0
    throw("pigpio initialisation failed.")
else
    @info "pigpio initialised okay."
end

channel = 0
baud = 32*10^3
spiFlags = 0b0000000000000000000000
spi = PiGPIOC.spiOpen(channel, baud, spiFlags) # channel, rate, flag

# init
send_16(spi, [0b0000_1100, 0b0000_0001]) # shut down off (normal)
send_16(spi, [0b0000_1111, 0b0000_0000]) # display test off (normal)
send_16(spi, [0b0000_1011, 0b0000_0111]) # no scan limit
send_16(spi, [0b0000_1001, 0b0000_0000]) # decode on for digits 0-7

for i in 1:10
sleep(2)
clear_all(spi)

# H
send_16(spi, [0b0000_1000, 0b0011_0111])
# E
send_16(spi, [0b0000_0111, 0b0100_1111])
# L
send_16(spi, [0b0000_0110, 0b0000_1110])
# L
send_16(spi, [0b0000_0101, 0b0000_1110])
# O
send_16(spi, [0b0000_0100, 0b0111_1110])

sleep(2)
clear_all(spi)

# J
send_16(spi, [0b0000_1000, 0b0011_1100])
# U
send_16(spi, [0b0000_0111, 0b0011_1110])
# L
send_16(spi, [0b0000_0110, 0b0000_1110])
# I
send_16(spi, [0b0000_0101, 0b0011_0000])
# A
send_16(spi, [0b0000_0100, 0b0111_0111])

end

println("Stop.")


text:
"""
# ArrayStream Construction

## Array -> ArrayStream

# Arguments

* `s::IO`: should be any subtype of IO for reading/writing data
* `streamindices::Union{AbstractRange,Int,Vector}`: linear indices of the `s`.
    This indicates which positions of `s` are to be read. LinearIndices (derived
    from the `S` trait, or the size) are mapped to `streamindices` to ensure
    appropriate indexing.
* `owsntream::Bool`: whether a given instance of `ArrayStream is the primary accesor to `s`
* `needswap::Bool`: indicates whether byte swapping is necessary at read time.

"""
struct ArrayStream{S,T,N,I<:Union{AbstractRange,Int,Vector}} <: AbstractArrayStream{S,T,N,I}
    s::IO
    streamindices::I
    ownstream::Bool
    needswap::Bool
end

"""
    read(ArrayStream)
    read!(ArrayStream, sink)
    read(AxisStream)
    read!(AxisStream, sink)
    read(MetaAxisStream)
    read!(MetaAxisStream, sink)

"""
# TODO: would it be better to use AbstractArrayStream here?
const AxisStream{S,T,N,I,Ax} = AxisArray{T,N,ArrayStream{S,T,N,I},Ax}
const MetaStream{S,T,N,I} = ImageMeta{T,N,ArrayStream{S,T,N,I}}
const MetaAxisStream{S,T,N,I,Ax} = ImageMeta{T,N,AxisStream{S,T,N,I,Ax}}
const MetaAxis{T,N,D,Ax} = ImageMeta{T,N,AxisArray{T,N,D,Ax}}
const AbstractMetaStream{S,T,N,I} = Union{MetaStream{S,T,N,I},MetaAxisStream{S,T,N,I,Ax}} where {Ax}
const AbstractStreamContainer = Union{AbstractMetaStream, AxisStream}

# TODO: maybe use Base.unsafe_wrap here
ArrayStream{S,T,N}(io::IO, streamindices::I, ownstream::Bool=true, needswap::Bool=false) where {S,T,N,I<:Union{AbstractRange,Int,Vector}} =
    ArrayStream{S,T,N,I}(io, streamindices, ownstream, needswap)

function ArrayStream(A::AbstractArray{T,N}, needswap::Bool=false) where {T,N}
    io = IOBuffer()
    write(io, A)
    seek(io, 0)
    ArrayStream{Tuple{size(A)...},T,N}(io, range(position(io), step=sizeof(T), length=length(A)), 0, true, needswap)
end

########
# Read #
########
function Base.read!(s::ArrayStream{S,T,N,AbstractRange}, sink::Array{T,N}) where {S,T,N}
    seek(s, firstindex(s))
    fixendian(read!(stream(s), sink), needswap(s))
end
Base.read(s::ArrayStream{S,T,N}) where {S,T,N} = read!(s, Array{T,N}(undef, size(s)))

# This is really slow
function Base.read!(s::ArrayStream{S,Tr,N,Vector{Int}}, sink::Array{Ts,N}) where {S,Tr<:BitTypes,Ts<:BitTypes,N}
    for i in 1:length(s)
        seek(s, streamindices(s)[i])
        sink[i] = read(stream(s), T)
    end
    return sink
end

# this needs testing but may be better
function Base.read!(s::ArrayStream{S,Tr,N,Vector{Int}}, sink::Array{Ts,N}) where {S,Tr<:BitTypes,Ts<:BitTypes,N}
    for i in 1:length(s)
        seek(s, streamindices(s)[i])
        sink[i] = read(stream(s), T)
    end
    return sink
end

Base.read!(s::ArrayStream{S,Tr,N}, sink::Array{Ts,N}) where {S,Tr<:BitTypes,Ts<:BitTypes,N} =
    fill!(sink, reinterpret(Ts, fixendian(sink, read(s))))

#########
# Write #
#########
Base.write(s::S, sink::AbstractArray) where {S<:AbstractStreamContainer} = write(getfield(s, :data), sink)
Base.write(s::ArrayStream, sink::AxisArray) = write(s, getfield(sink, :data))
Base.write(s::ArrayStream, sink::ImageMeta) = write(s, getfield(sink, :data))

function Base.write(s::ArrayStream{S,T,N,AbstractRange}, sink::Array{T,N}) where {S,T,N}
    seek(s, firstindex(s))
    write(stream(s), fixendian(sink, needswap(s)))
end

function Base.read!(s::ArrayStream{S,T,N,Vector{Int}}, sink::Array{T,N}) where {S,T,N}
    for i in 1:length(s)
        seek(s, streamindices(s)[i])
        write(stream(s), fixendian(sink[i]), needswap(s))
    end
end

Base.write(s::ArrayStream{S,<:BitTypes,N}, sink::Array{<:BitTypes,N}) where {S,N} =
    write(s, reinterpret(Tr, s))


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Data: 06/04/2014

module Util

include("format_string.jl")

export outputGFormatString

end

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
export Gene

mutable struct Gene
    inno_nb::Integer
    origin::Float64
    destination::Float64
    weight::Float64
    enabled::Bool
end

function Gene(origin::Float64, destination::Float64, inno::Integer)
    w = rand_weight()
    Gene(inno, origin, destination, w, true)
end

function Gene(g::Gene)
    deepcopy(g)
end

function Gene(d::Dict{String,Any})
    Gene(
        d["inno_nb"],
        d["origin"],
        d["destination"],
        d["weight"],
        d["enabled"],
    )
end


text:
const domainXY = range(-maxR, stop = maxR, length = 333)/w₀
const domainR  = range(0., stop = maxR, length = 800)/w₀
const time = range(- n * τ, stop = n * τ, length = 1200)/T₀

# useful function definitions

function w_inst(x, y, z, t, type)
    @unpack ε₀, c = fundamental_constants(type)
    return (ε₀/2)*(E([x,y,z], t, type)⋅E([x,y,z], t, type) +
               c^2*B([x,y,z], t, type)⋅B([x,y,z], t, type))
end

# dot takes the conjugate of the first vector automatically!
function w_cycle(x, y, z, type)
    @unpack ε₀, c = fundamental_constants(type)
    return (ε₀/4)*real(E([x, y, z], type)⋅E([x, y, z], type) +
                 (c^2)*B([x, y, z], type)⋅B([x, y, z], type))
end

function j_inst(x, y, z, t, type)
    ε₀ = fundamental_constants(type, :ε₀)
    return ε₀*([x,y,z]×(E([x,y,z], t, type)×B([x,y,z], t, type)))
end

function jz_cycle(x, y, z, type)
    ε₀ = fundamental_constants(type, :ε₀)

    return (ε₀/4)*real((x*(
                            (E([x, y, z], type)×conj(B([x, y, z], type)))[2] +
                            (conj(E([x, y, z], type))×B([x, y, z], type))[2]
                            )
                            -
                        y*(
                            (E([x, y, z], type)×conj(B([x, y, z], type)))[1] +
                            (conj(E([x, y, z], type))×B([x, y, z], type))[1]
                            ))
                            )
end

function Mz(x, y, z, t, type) # normed to the q*E₀*w₀
    q = type.constants.q
    E₀ = type.derived.E₀
    w₀ = type.w₀
    sign(q)*(x*E([x, y, z], t, type)[2] - y*E([x, y, z], t, type)[1])/(E₀*w₀)
end

E_ρAvg(ρ, z, t, type) = quadgk(ϕ -> E([ρ*cos(ϕ), ρ*sin(ϕ), z], t, type)/(2*π), 0, 2*π, rtol = 1e-8, atol = 1e-10)[1]

wi_ρAvg(ρ, z, t, type) = quadgk(ϕ -> w_inst(ρ*cos(ϕ), ρ*sin(ϕ), z, t, type)/(2*π), 0, 2*π, rtol = 1e-8)[1]
wc_ρAvg(ρ, z, type) = quadgk(ϕ -> w_cycle(ρ*cos(ϕ), ρ*sin(ϕ), z, type)/(2*π), 0, 2*π, rtol = 1e-8)[1]

ji_ρAvg(ρ, z, t, type) = quadgk(ϕ -> j_inst(ρ*cos(ϕ), ρ*sin(ϕ), z, t, type)/(2*π), 0, 2*π, rtol = 1e-8, atol = 1e-10)[1]
jc_ρAvg(ρ, z, type) = quadgk(ϕ -> jz_cycle(ρ*cos(ϕ), ρ*sin(ϕ), z, type)/(2*π), 0, 2*π, rtol = 1e-8, atol = 1e-10)[1]

# end of definitions

common_settings = (fontfamily = "computer modern", dpi = 360, label = false)
plot_settings   = (color_palette = :Set1_4, seriescolor = 4, linewidth = 2.5)

for type in lasers
    pol = polarization(type)

    folder_setup("$(pol)_p$(type.p)_m$(type.m)") # change folder for each simulation
    folder_setup("fields") # separate fields from particles

    # w₀ = type.w₀
    p = type.p
    m = type.m
    ξy = type.polarization.ξy
    # T₀ = type.derived.T₀
    roots = rootsLUT[(p = p, m = abs(m))]

    file_ext = "$(pol)_p$(p)_m$(m).$(fmt)"

    # 1. Heatmaps

    w_00 = heatmap(domainXY, domainXY, (x,y) -> w_inst(w₀*x, w₀*y, 0, 0, type);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :linear_kryw_5_100_c67_n256,
                clims = x -> (0, maximum(x)),
                aspect_ratio = 1, common_settings...)
    add_info!(w_00, type, true)
    savefig(w_00, "xy00Energy_$file_ext")
    println("Energy density for time = 0, z = 0, in the xy-plane")

    w_0T₀ = heatmap(domainXY, domainXY, (x,y) -> w_cycle(w₀*x, w₀*y, 0, type);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :linear_kryw_5_100_c67_n256,
                clims = x -> (0, maximum(x)),
                aspect_ratio = 1, common_settings...)
    add_info!(w_0T₀, type, true)
    savefig(w_0T₀, "xy0T₀Energy_$file_ext")
    println("Energy density cycle averaged, z = 0, in the xy-plane")

    Ex_00 = heatmap(domainXY, domainXY, (x,y) -> E([w₀*x, w₀*y, 0], 0, type)[1];
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(Ex_00, type, true)
    savefig(Ex_00, "xy00Ex_$file_ext")
    println("Electric field on the x-axis for time = 0, z = 0, in the xy-plane")

    Ez_00 = heatmap(domainXY, domainXY, (x,y) -> E([w₀*x, w₀*y, 0], 0, type)[3];
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(Ez_00, type, true)
    savefig(Ez_00, "xy00Ez_$file_ext")
    println("Electric field on the z-axis for time = 0, z = 0, in the xy-plane")

    jz_00 = heatmap(domainXY, domainXY, (x,y) -> j_inst(w₀*x, w₀*y, 0, 0, type)[1];
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(jz_00, type, true)
    savefig(jz_00, "xy00jz_$file_ext")
    println("z-axis field angular momentum density for time = 0, z = 0, in the xy-plane")

    jz_0T₀ = heatmap(domainXY, domainXY, (x,y) -> jz_cycle(w₀*x, w₀*y, 0, type);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(jz_0T₀, type, true)
    savefig(jz_0T₀, "xy0T₀jz_$file_ext")
    println("z-axis field angular momentum density cycle averaged, z = 0, in the xy-plane")

    jz_w = heatmap(domainXY, domainXY, (x,y) -> ω*jz_cycle(w₀*x, w₀*y, 0, type)/w_cycle(w₀*x, w₀*y, 0, type);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :lajolla,
                aspect_ratio = 1, common_settings...)
    add_info!(jz_w, type, true)
    savefig(jz_w, "xy0T₀jz_w_$file_ext")
    println("z-axis field angular momentum density divided by energy density, cycle averaged, z = 0, in the xy-plane")

    Lz′ = heatmap(domainXY, domainXY, (x,y) -> Mz(w₀*x, w₀*y, 0, 0, type);
                xlabel = L"x/w_0", ylabel = L"y/w_0",
                seriescolor = :vik,
                clims = x -> maximum(abs.(extrema(x))).*(-1,1),
                aspect_ratio = 1, common_settings...)
    add_info!(Lz′, type, true)
    savefig(Lz′, "xy00Lz_$file_ext")
    println("Instantaneous reduced torque for time t = 0, z = 0, in the xy-plane")

    @time begin
    Lz′_avg = heatmap(domainXY, domainXY, (x,y) -> quadgk(t -> Mz(w₀*x, w₀*y, 0, t, type), -Inf, Inf, atol = 1e-2)[1];
    xlabel = L"x/w_0", ylabel = L"y/w_0",
    seriescolor = :vik,
    clims = x -> maximum(abs.(extrema(x))).*(-1,1),
    aspect_ratio = 1, common_settings...)
    add_info!(Lz′_avg, type, true)
    savefig(Lz′_avg, "xy0∞Lz_$file_ext")
    println("Pulse averaged reduced torque for, z = 0, in the xy-plane")
    end
    # 2. Azimuthally averaged radial plots

    jz_ρAvg = plot(domainR, r -> jc_ρAvg(w₀*r, 0, type);
                label = false)
    jz_ρAvg = add_rect!(jz_ρAvg, domainR, r -> jc_ρAvg(w₀*r, 0, type), roots;
                xlabel = L"\rho/w_0", ylabel = L"\overline{j_z^{xy}}(\rho, t = 0)\quad\mathrm{[a.u.]}",
                common_settings..., plot_settings...)
    add_info!(jz_ρAvg, type)
    savefig(jz_ρAvg, "jz_ρAvg_$file_ext")
    println("Azimuthally averaged angular momentum on the z-axis for time = 0, z = 0, in the xy")

    w_rAvg = plot(domainR, r -> wc_ρAvg(w₀*r, 0, type);
                label = false)
    w_rAvg = add_rect!(w_rAvg, domainR, r -> wc_ρAvg(w₀*r, 0, type), roots;
        xlabel = L"\rho/w_0", ylabel = L"\mathrm{Intensity}(\rho)\quad\mathrm{[a.u.]}",
        common_settings..., plot_settings...)
    add_info!(w_rAvg, type)
    savefig(w_rAvg, "w_ρAvg_$file_ext")
    println("Azimuthally averaged energy density for time = 0, z = 0, in the xy")

    # 3. Time axis plots.

    Ex_t = plot(time, t -> E([w₀/√2/2, w₀/√2/2, 0], t*T₀, type)[1];
         xlabel = L"t/T_0", ylabel = L"E_x({w_0}/{\sqrt{8}}, {w_0}/{\sqrt{8}}, 0, t)\quad\mathrm{[a.u.]}",
         common_settings..., plot_settings...)
    add_info!(Ex_t, type)
    savefig(Ex_t,"Ex_t_$file_ext")
    println("Electric field on the x-axis in time for x = w₀/√2/2, y = w₀/√2/2, z = 0")

    # 3.1 Animations

    animLz = Animation()
    for (i,t) in enumerate(range(-n*T₀, n*T₀, length = 299)) # +1 = 300 / 30fps = 10s
        Lz′_anim = heatmap(domainXY, domainXY, (x,y) -> Mz(w₀*x, w₀*y, 0,  t, type);
                    xlabel = L"x/w_0", ylabel = L"y/w_0",
                    seriescolor = :vik,
                    clims = x -> (-0.5,0.5), # maximum(abs.(extrema(x))).*(-1,1),
                    cbar = :none,
                    aspect_ratio = 1,
                    fontfamily = "computer modern",
                    dpi = 120,
                    label = false)
        add_info!(Lz′_anim, type, true)
        frame(animLz)
        pritnln("Animatin ... Frame $i/300.")
    end
    gif(animLz, "Lz′.gif", fps = 30)
    println("Instantaneous reduced torque for time t = 0, z = 0, in the xy-plane")

    cd("../..") # return to output_* base folder
    break
end

max_stars_repo_name: toma-vlad/electrodynamic-simulations

max_issues_repo_name: toma-vlad/electrodynamic-simulations

max_forks_repo_name: toma-vlad/electrodynamic-simulations


text:
## constructors.jl
## (c) 2014--2020 David A. van Leeuwen

## Constructors related to the types in namedarraytypes.jl

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

letter(i) = string(Char((64+i) % 256))

## helpers for constructing names dictionaries
defaultnames(dim::Integer) = map(string, 1:dim)
function defaultnamesdict(names::AbstractVector)
    dict = OrderedDict(zip(names, 1:length(names)))
    length(dict) == length(names) || error("Cannot have duplicated names for indices")
    return dict
end
defaultnamesdict(dim::Integer) = defaultnamesdict(defaultnames(dim))
defaultnamesdict(dims::Tuple) = map(defaultnamesdict, dims) # ::NTuple{length(dims), OrderedDict}

defaultdimname(dim::Integer) = Symbol(letter(dim))
defaultdimnames(ndim::Integer) = ntuple(defaultdimname, ndim)
defaultdimnames(a::AbstractArray) = defaultdimnames(ndims(a))

## disambiguation (Argh...)
function NamedArray(a::AbstractArray{T,N},
                    names::Tuple{},
                    dimnames::NTuple{N, Any}) where {T,N}
    NamedArray{T,N,typeof(a),Tuple{}}(a, (), ())
end

NamedArray(a::AbstractArray{T,N}, names::Tuple{}) where {T,N} = NamedArray{T,N,typeof(a),Tuple{}}(a, (), ())
NamedArray(a::AbstractArray{T,0}, ::Tuple{}, ::Tuple{}) where T = NamedArray{T,0,typeof(a),Tuple{}}(a, (), ())

## Basic constructor: array, tuple of dicts, tuple of dimnames
function NamedArray(array::AbstractArray{T,N},
                    names::NTuple{N,OrderedDict},
                    dimnames::NTuple{N}=defaultdimnames(array)) where {T,N}
    ## inner constructor
    NamedArray{T, N, typeof(array), typeof(names)}(array, names, dimnames)
end

## constructor with array, names and dimnames (dict is created from names)
function NamedArray(array::AbstractArray{T,N},
                    names::NTuple{N,AbstractVector}=tuple((defaultnames(d) for d in size(array))...),
                    dimnames::NTuple{N, Any}=defaultdimnames(array)) where {T,N}
    dicts = defaultnamesdict(names)
    NamedArray{T, N, typeof(array), typeof(dicts)}(array, dicts, dimnames)
end

## Deprecated: use tuples, as above
## vectors instead of tuples, with defaults (incl. no names or dimnames at all)
function NamedArray(array::AbstractArray{T,N},
                    names::AbstractVector{VT},
                    dimnames::AbstractVector=[defaultdimname(i) for i in 1:ndims(array)]) where {T,N,VT}
    length(names) == length(dimnames) == N || error("Dimension mismatch")
    if VT == Union{} ## edge case, array == Array{}()
        dicts = ()
    elseif VT <: OrderedDict
        dicts = tuple(names...)
    else
        dicts = defaultnamesdict(tuple(names...))::NTuple{N, OrderedDict{eltype(VT),Int}}
    end
    NamedArray{T, N, typeof(array), typeof(dicts)}(array, dicts, tuple(dimnames...))
end

## I can't get this working
## @Base.deprecate NamedArray(array::AbstractArray{T,N}, names::AbstractVector{<:AbstractVector}, dimnames::AbstractVector) where {T,N} NamedArray(array::AbstractArray{T,N}, names::NTuple{N}, dimnames::NTuple{N}) where {T,N}

## special case for 1-dim array to circumvent Julia tuple-comma-oddity, #86
NamedArray(array::AbstractVector{T},
           names::AbstractVector{VT}=defaultnames(length(array)),
           dimname=defaultdimname(1)) where {T,VT} =
    NamedArray(array, (names,), (dimname,))

function NamedArray(array::AbstractArray{T,N}, names::NamedTuple) where {T, N}
    length(names) == N || error("Dimension mismatch")
    return NamedArray(array, values(names), keys(names))
end

## Type and dimensions
"""
`NamedArray(T::Type, dims::Int...)` creates an uninitialized array with default names
for the dimensions (`:A`, `:B`, ...) and indices (`"1"`, `"2"`, ...).
"""
NamedArray(T::DataType, dims::Int...) = NamedArray(Array{T}(undef, dims...))

NamedArray{T}(n...) where {T} = NamedArray(Array{T}(undef, n...))


text:
# Indexers
# ========
#
# Mapping from names to objects.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

module Indexers

export
    Indexer, names!, rename!, rename

"""
The Indexer type

The type is parametric and contains a dictionary mapping
names (Symbols) to objects, as well as a vector of names
(Symbols) that enables indexers to have order (i.e.
The user should be able to push and pop and insert with indexes).
"""
immutable Indexer{T}
    lookup::Dict{Symbol, T}
    names::Vector{Symbol}
end

"""
    Indexer{T}(::Type{T})

Construct an empty Indexer mapping Symbols to a given Type T.

# Examples
```julia
Indexer(Int64)
Indexer(Int8)
Indexer(UInt32)
```
"""
function Indexer{T}(::Type{T})
    return Indexer{T}(Dict{Symbol, T}(), Vector{Symbol}())
end

"""
    Indexer{T}(names::Vector{Symbol}, vals::Vector{T})

Create an index that associates a Vector of Symbols with a vector of values.

This constructor constructs the index such that names[1] -> vals[1],
names[2] -> vals[2] ... and so on.

# Examples
```julia
Indexer([:Hello, :Bye], [1, 2])
Indexer([:Hello, :Bye], [1:10, 11:20])
```
"""
function Indexer{T}(names::Vector{Symbol}, vals::Vector{T})
    u = make_unique(names)
    lookup = Dict{Symbol, T}(zip(u, vals))
    return Indexer{T}(lookup, u)
end

"""
    Indexer{S <: AbstractString, T}(names::Vector{S}, vals::Vector{T})

Create an index that associates names in a Vector of Strings with a vector of
values.

This constructor constructs the index such that names[1] -> vals[1],
names[2] -> vals[2] ... and so on.

# Examples
```julia
Indexer(["Hi", "Bye", "Sup"], [1, 2, 3])
Indexer(["Hi", "Bye", "Sup"], [1:5, 6:10, 11:15])
```
"""
function Indexer{S <: AbstractString, T}(names::Vector{S}, vals::Vector{T})
    Indexer(convert(Vector{Symbol}, names), vals)
end

"""
    Indexer(names::Vector{Symbol}, ::Type{T})

Construct an Indexer mapping a Vector of Symbols (names), to a set of
values of type T.

Provide this constructor a vector of Symbols as names, and a type that is a
subtype of Integer, and it will construct an index such that
names[1] -> 1, names[2] -> 2, ..., names[length(names)] -> length(names).

# Examples
```julia
Indexer([:First, :Second, :Third], Int64)
Indexer([:First, :Second, :Third], UInt32)
```
"""
function Indexer{T <: Integer}(names::Vector{Symbol}, ::Type{T})
    return Indexer(names, collect(T(1):T(length(names))))
end

"""
    Indexer{S <: AbstractString, T <: Integer}(names::Vector{S}, ::Type{T})

Construct an Indexer mapping a Vector of Strings (names), to a set of
values of type T <: Integer.

Provide this constructor a vector of Strings as names, and a type that is a
subtype of Integer, and it will construct an index such that
names[1] -> 1, names[2] -> 2, ..., names[length(names)] -> length(names).

# Examples
```julia
Indexer(["First", "Second", "Third"], Int64)
Indexer(["First", "Second", "Third"], UInt32)
```
"""
function Indexer{S <: AbstractString, T <: Integer}(names::Vector{S}, ::Type{T})
    return Indexer(convert(Vector{Symbol}, names), T)
end

"""
    Indexer(names::Vector{Symbol})

Construct an Indexer mapping a Vector of Symbols (names), to a set of
values of type Int, where Int is the default integer encoding on your machine
(usually Int32 or Int64).

Provide this constructor a vector of Symbols as names, and it will construct
an index such that names[1] -> 1, names[2] -> 2, ...,
names[length(names)] -> length(names).

# Examples
```julia
Indexer([:First, :Second, :Third])
```
"""
function Indexer(names::Vector{Symbol})
    return Indexer(names, Int)
end

"""
    Indexer{S <: AbstractString}(names::Vector{S})

Construct an Indexer mapping a Vector of Strings (names), to a set of
values of type Int, where Int is the default integer encoding on you machine
(usually Int32 or Int64).

Provide this constructor a vector of Strings as names, and it will construct
an index such that names[1] -> 1, names[2] -> 2, ...,
names[length(names)] -> length(names).

# Examples
```julia
Indexer(["First", "Second", "Third"])
```
"""
function Indexer{S <: AbstractString}(names::Vector{S})
    return Indexer(convert(Vector{Symbol}, names), Int)
end

"""
    make_unique(names::Vector{Symbol})

Make sure all vectors in a symbol are unique.

# Examples
```julia
make_unique([:First, :Second, :Third])
make_unique([:First, :First, :First])
make_unique([:First, :Third, :Third])
```
"""
function make_unique(names::Vector{Symbol})
    seen = Set{Symbol}()
    names = copy(names)
    dups = Int[]
    for (i, name) in enumerate(names)
        in(name, seen) ? push!(dups, i) : push!(seen, name)
    end
    for i in dups
        nm = names[i]
        k = 1
        while true
            newnm = Symbol("$(nm)_$k")
            if !in(newnm, seen)
                warn("Trying to build an index with duplicate names. Changing $(names[i]) to $(newnm)")
                names[i] = newnm
                push!(seen, newnm)
                break
            end
            k += 1
        end
    end
    return names
end

Base.length(x::Indexer) = length(x.names)
Base.names(x::Indexer) = copy(x.names)
Base.copy(x::Indexer) = Indexer(copy(x.lookup), copy(x.names))
Base.:(==)(x::Indexer, y::Indexer) = (x.lookup == y.lookup) && (x.names == y.names)
Base.haskey(x::Indexer, key::Symbol) = haskey(x.lookup, key)
Base.haskey(x::Indexer, key::AbstractString) = haskey(x, convert(Symbol, key))
Base.haskey(x::Indexer, key::Integer) = 1 <= key <= length(x.names)
Base.keys(x::Indexer) = names(x)

"""
    names!(x::Indexer, names::Vector{Symbol})

Completely replace the names in an Indexer with a new Vector of Symbols.

# Examples
```julia
names!(my_indexer, [:First, :NewName, :Third])
```
"""
function names!(x::Indexer, names::Vector{Symbol})
    if length(names) != length(x)
        throw(ArgumentError("Length of new names doesn't match length of Index."))
    end
    names = make_unique(names)
    if length(x.names) != length(x.lookup)
        throw(ArgumentError("Index corrupted by duplicate symbols in names."))
    end
    rename!(x, x.names, names)
    return x
end

"""
    names!{S <: AbstractString}(x::Indexer, names::Vector{S})

Completely replace the names in an Indexer with a new Vector of Strings.

# Examples
```julia
names!(my_indexer, ["First", "NewName", "Third"])
```
"""
function names!{S <: AbstractString}(x::Indexer, names::Vector{S})
    names!(x, convert(Vector{Symbol}, names))
end

function rename!(x::Indexer, names)
    for (from, to) in names
        if !haskey(x, from)
            throw(ArgumentError("Cannot rename $from to $to, $from does not exist in the Indexer."))
        end
        if haskey(x, to)
            throw(ArgumentError("Cannot rename $from to $to, $to already exists in the Indexer."))
        end
        # Change the name array
        nameidx = findfirst(x.names, from)
        x.names[nameidx] = to
        # Change the dictionary
        x.lookup[to] = val = pop!(x.lookup, from)
    end
    return x
end

"""
    rename!(x::Indexer, from::Vector{Symbol}, to::Vector{Symbol})

Rename a set of names in an Indexer, to a new set of names.

This will rename such that the name in from[1] will be renamed to to[1],
from[2] will be renamed to to[2] and so on.

# Examples
```julia
rename!(my_indexer, [:First, :Third], [:NewFirst, :NewThird])
```
"""
function rename!(x::Indexer, from::Vector{Symbol}, to::Vector{Symbol})
    return rename!(x, zip(from, to))
end

"""
    rename!{S <: AbstractString}(x::Indexer, from::Vector{S}, to::Vector{S})

Rename a set of names in an Indexer, to a new set of names.

This will rename such that the name in from[1] will be renamed to to[1],
from[2] will be renamed to to[2] and so on.

# Examples
```julia
rename!(my_indexer, ["First", "Third"], ["NewFirst", "NewThird"])
```
"""
function rename!{S <: AbstractString}(x::Indexer, from::Vector{S}, to::Vector{S})
    return rename!(x, convert(Vector{Symbol}, from), convert(Vector{Symbol}, to))
end

"""
    rename!(x::Indexer, from::Symbol, to::Symbol)

Rename a single name (from) in Indexer x, to a new name (to).

# Examples
```julia
rename!(my_indexer, :First, :NewFirst)
```
"""
rename!(x::Indexer, from::Symbol, to::Symbol) = rename!(x, ((from, to),))

"""
    rename!{S <: AbstractString}(x::Indexer, from::S, to::S)

Rename a single name (from) in Indexer x, to a new name (to).

# Examples
```julia
rename!(my_indexer, "First", "NewFirst")
```
"""
function rename!{S <: AbstractString}(x::Indexer, from::S, to::S)
    return rename!(x, ((convert(Symbol, from), convert(Symbol, to)),))
end

"""
    rename!(f::Function, x::Indexer)

Rename a indexer according to some function, f, that accepts one argument.

Each name in the Indexer x, will be input to, f, and the output will be
what that name will be renamed to.
Therefore the function, f, will need to return Symbols or Strings.
"""
rename!(f::Function, x::Indexer) = rename!(x, [(x,f(x)) for x in x.names])
rename(x::Indexer, args...) = rename!(copy(x), args...)

# Indexing into an index mapping single names to single numbers

## Indexing with single name
Base.getindex(x::Indexer, key::Symbol) = x.lookup[key]
Base.getindex(x::Indexer, key::AbstractString) = x.lookup[convert(Symbol, key)]

## Indexing with multiple names.
Base.getindex(x::Indexer, keys::Vector{Symbol}) = [x.lookup[key] for key in keys]
function Base.getindex{S <: AbstractString}(x::Indexer, keys::Vector{S})
    return [x.lookup[convert(Symbol, key)] for key in keys]
end

## Indexing with a single integer
function Base.getindex{I <: Integer}(x::Indexer, key::I)
    n = x.names[key]
    return x[n]
end

## Indexing with a range of integers
function Base.getindex{I <: Integer}(x::Indexer, keys::AbstractVector{I})
    ns = x.names[keys]
    return x[ns]
end

## Indexing with a vector of bools.
function Base.getindex(x::Indexer, keys::Vector{Bool})
    ns = x.names[keys]
    return x[ns]
end

end


text:
struct Met{T<:Real}
    id::String
    S::Vector{T}
    rxns::Vector
    b::T

    function Met{T}(id::String;
            S = T[],
            rxns = [],
            b = zero(T)
        ) where {T<:Real}

        length(S) != length(rxns) && error("'S' and 'rxns' must have the same length")
        return new{T}(id, S, rxns, b)
    end
end

Met(id::String; kwargs...) = Met{Float64}(id; kwargs...)


text:
using DataStructures

const score = Dict(
  ')' => 3,
  ']' => 57,
  '}' => 1197,
  '>' => 25137
)

const score2 = Dict(
  ')' => 1,
  ']' => 2,
  '}' => 3,
  '>' => 4
)

const closing = Dict(
  '(' => ')',
  '[' => ']',
  '{' => '}',
  '<' => '>'
)

const openChar  = Set(['(', '[', '{', '<'])
const closeChar = Set([')', ']', '}', '>'])

# Calculates the score for problem 2
function calcScore(s::Array{Char}) :: Int
  res = 0
  map( x -> res = 5*res + score2[x], s)
  return res
end

# Evaluates expression, and returns the following code:
#   1  -> Expression is valid and complete
#   0  -> Expression is valid but incomplete
#   -1 -> Expression is invalid
# If invalid, returns invalid character.
# If incomplete, returns array of missing characters.
function eval(expr::String) :: Tuple{Int,Vector{Char}}
  s = Stack{Char}()
  for c in expr
    in(c,closeChar) && ( isempty(s) || (c != closing[pop!(s)])) && return (-1,[c]) # Closing char
    in(c, openChar) && push!(s,c) # Opening char, put in stack
  end
  !isempty(s) && return (0,map(c -> closing[c], s))
  return (1,[])
end

function readInput() :: Vector{String}
  f = open("./inputs/day10.txt")
  data :: Vector{String} = []

  for line in eachline(f)
    push!(data, line)
  end
  return data
end

function problem1(lines::Vector{String}) :: Int
  return sum(map(x -> score[x[2][1]] , filter(x -> x[1] == -1,map(line -> eval(line), lines))))
end

function problem2(lines::Vector{String}) :: Int
  scores = sort(map(x -> calcScore(x[2]) , filter(x -> x[1] == 0,map(line -> eval(line), lines))))
  return scores[length(scores) ÷ 2 + 1]
end

lines = readInput()

println("Problem 1:")
res = problem1(lines)
println(res)

println("Problem 2:")
res = problem2(lines)
println(res)


text:
using SMACOF
using Plots
using Random
gr()
ENV["GKSwstype"] = "100"
Random.seed!(2021)

include("../example_template.jl")

gifdir = "examples/wda_smacof/"
filename = "noiseless"

generate_gif((X, Δ) -> wda_smacof(Δ,
                    η=0.9,
                    Xinit=rand(size(X)...),
                    verbose=true,
                    anchors=X,
                    ε=1e-4,
                    DA_itmax=30,
                    return_history=true)[2],
            gifdir, filename, n=500)


text:
using StatsBase

abstract StepData

immutable Outputter
	outputsFunction::Function
	outputStringHeader::AbstractString
	numOutputsFromOutputsFunction::Int64
end

immutable ResultFromOO
	resultString::AbstractString
	resultLine::Vector{Float64}
end

immutable OutputOpts
	ooString::AbstractString
	outputter::Outputter
	expIndices
	function OutputOpts(outputter::Outputter,expIndices; average::Bool=false, ooString::AbstractString="average = $average")
		new(ooString, outputter,expIndices)
	end
end

immutable Opts
	init::Matrix{Float64}
	stepSizePower::Int64
	maxG::Int64
	outputLevel::Int64
	stepOutputLevel::Int64
	optsString::AbstractString
	function Opts(init::Matrix{Float64}; stepSizePower::Int64=1, maxG::Int64=typemax(Int64)-35329, outputLevel::Int64=1, optsString::AbstractString="stepSizePower = $stepSizePower, maxG = $maxG", stepOutputLevel::Int64=0)
		outputLevel>2 && println("outputLevel is $outputLevel")
		(maxG>typemax(Int64)-35329 || maxG<0) && error("maxG is $maxG, and is out of range")
		new(init, stepSizePower, maxG, outputLevel, stepOutputLevel, optsString)
	end
end

immutable Problem
	L2reg::Bool
	name::AbstractString
	lossFunctionString::AbstractString
	numTP::Int64 # Don't need for EGR
	getNextSampleFunction::Task
	getSampleFunctionAt::Function # Don't need for EGR
	function Problem(a,b,c,d,e,f)
		new(a,b,c,d,e,f)
	end
end

function alg(problem::Problem, opts::Opts, sd::StepData, oo::OutputOpts, writeFunction::Function)

	srand(1)

	opts.outputLevel>0 && println("Function to minimize: $(problem.name) with loss function $(problem.lossFunctionString), L2reg = $(problem.L2reg)")
	opts.outputLevel>0 && println("Step: $(sd.stepString) with stepSizePower = $(opts.stepSizePower)")
	opts.outputLevel>0 && println("Opts: $(opts.optsString)")
	opts.outputLevel>0 && println("OutputOpts: $(oo.ooString)")

	maxG = min(opts.maxG, oo.expIndices[end]);

	results_k = Int64[]
	results_gnum = Int64[]
	results_fromOutputsFunction = zeros(0,oo.outputter.numOutputsFromOutputsFunction)

	kOutputs=1
	k=0
	gnum=0
	x=opts.init

	xSum=copy(x)
	xToTest = zeros(size(x))

	opts.outputLevel>1 && println("        k     gnum"*oo.outputter.outputStringHeader)

	while true

		if gnum >= oo.expIndices[kOutputs]
			xToTest =x
			fromOutputsFunction::ResultFromOO = oo.outputter.outputsFunction(xToTest)
			if opts.outputLevel>1
				@printf("%2.i %8.i %8.i ", kOutputs, k, gnum)
				println(fromOutputsFunction.resultString)
			end
			writeFunction(problem, sd, opts, k, gnum, oo.expIndices[kOutputs], fromOutputsFunction)
			if isnan(fromOutputsFunction.resultLine[end])
				return ("NaN found", results_k, results_gnum,results_fromOutputsFunction,xToTest)
			end
			push!(results_k,k)
			push!(results_gnum,gnum)
			results_fromOutputsFunction = [results_fromOutputsFunction ; fromOutputsFunction.resultLine']
			kOutputs += 1
		end

		gnum >= maxG && break

		(g, gnum) = sd.getStep(x, k, gnum, sd, problem; outputLevel = opts.stepOutputLevel)

		opts.outputLevel>2 && println("norm(g) before step = $(norm(g))")
		opts.outputLevel>2 && println("norm(x) before step = $(norm(x))")

		x-=(2.0^opts.stepSizePower)*g

		opts.outputLevel>2 && println("norm(x) after step = $(norm(x))")

		xSum += x;

		k+=1

	end

	# exit()

	opts.outputLevel>0 && println("Finished alg: $(sd.stepString)")

	## Write answer!
	("Finished nicely", results_k, results_gnum, results_fromOutputsFunction, xToTest)
end


text:
module RayTracer

export  Point , Vec , Vec1 , Vec2 , Vec3
export Color
export Canvas
export Float

using Parameters , StaticArrays

const Float = Float64

include("Vec.jl")
include("Point.jl")
include("Color.jl")
include("Canvas.jl")
end


text:
# _updaterules
# Update the StaticRuleset in the SimData object with
# a (potentially) modifed version of the original rules tuple.
# This must be passed in to allow async updating from a live
# interface while the simulation runs.
function _updaterules(rules::Tuple, sd::AbstractSimData)
    newrs = ModelParameters.setparent(
        ruleset(sd),
        _modifyrules(_proc_setup(proc(sd), ModelParameters.stripparams(rules)), sd)
    )
    @set sd.ruleset = newrs
end

# _modifyrules
# Run `modifyrule` for each rule, recursively.
_modifyrules(rules::Tuple, simdata) =
    (modifyrule(rules[1], simdata), _modifyrules(tail(rules), simdata)...)
_modifyrules(rules::Tuple{}, simdata) = ()

# The default `modifyrule` returns the rule unchanged.
modifyrule(rule, simdata) = rule