text:
##############################
## by Qin Yu, UCL, Nov 2018
## using Julia 1.0.1
##############################

############################## Starting here:
using LinearAlgebra
using Statistics
using Random
using Plots
using MAT  # Windows only
using JLD  # macOS Friendly
using CSV
using Printf
using DataFrames
push!(LOAD_PATH, ".")
using SuperLearn
gr()

############################## The Boston Housing Dataset:
## Detail see 3.Linear_Regression_Boston_Housing.jl

############################## Import Boston:
#----- (for my Windows 10)
# Import Boston:
matdata = matread("boston.mat")["boston"]
# If MAT package is not available for your version of Julia,
# use JLD to prepare data on another machine:
# save("./boston.jld", "boston", matdata)

#----- (for my macOS)
# use JLD to load the transfored data:
matdata = load("boston.jld")["boston"]
# and have a look at all data of y:
# plot(matdata[:, 14])

############################## 5-fold Cross-validation to Find Best Parameter Pair
# Define a pool of parameters, we want to find the best pair:
𝛄 = [2.0^i for i in -40:-26]  # regularisation parameter γ, and
𝛔 = [2.0^i for i in collect(7:.5:13)]  # variance parameter σ for Gaussian kernel

# Prepare 5-folds data:
nrow, ncol = size(matdata)
nrow_test  = div(nrow, 3)
nrow_train = nrow - nrow_test
Random.seed!(1111)  # this number is again pure magic
permuted_matdata = matdata[randperm(nrow),:]
permuted_matdata_train = permuted_matdata[1:nrow_train,:]
permuted_matdata_test = permuted_matdata[nrow_train+1:end,:]

nrow_5fold  = div(nrow_train, 5)
𝑙 = nrow_train - nrow_5fold
𝑰 = Matrix{Float64}(I, 𝑙, 𝑙)  # 𝑰 = Id(l by l)

𝐾(𝒙ᵢ, 𝒙ⱼ, σ) = exp(-(norm(𝒙ᵢ - 𝒙ⱼ)^2)/(2 * σ^2))  # 𝐾(𝒙ᵢ, 𝒙ⱼ) = 𝑒^{∥𝒙ᵢ-𝒙ⱼ∥²/2σ²}
function fill_𝑲(𝑲, 𝑿, σ)
    nrow_𝑲, ncol_𝑲 = size(𝑲)
    for i in 1:nrow_𝑲
        for j in 1:ncol_𝑲
            𝑲[i, j] = 𝐾(𝑿[i,:], 𝑿[j,:], σ)  # 𝐾ᵢⱼ = 𝐾(𝒙ᵢ, 𝒙ⱼ)
        end
    end
end
get_𝜶(𝑲, γ, 𝑰, 𝒚) = inv(𝑲 + (γ * 𝑙 * 𝑰)) * 𝒚  # 𝛂* = (𝑲 + γ𝑙𝑰)⁻¹⋅𝒚
function get_ŷ(𝜶, x, σ, γ, 𝑿_train)  # 𝒙ₐ denoted by x
    𝐾ₐ(𝒙ᵢ, σ) = 𝐾(𝒙ᵢ, x', σ)
    𝑿_train_vectorised = [𝑿_train[i,:]' for i = 1:size(𝑿_train, 1)]
    ŷ = 𝜶' * 𝐾ₐ.(𝑿_train_vectorised, σ)  # ŷ = 𝜶ᵀ⋅(𝐾₁ₐ, 𝐾₂ₐ, ..., 𝐾ₗₐ)ᵀ
end

SSE = MSE = zeros(size(𝛔,1), size(𝛄,1))
# DO NOT RUN THE TRI-LOOP!!! because otherwise you need to wait.
# If you don't want to wait, call these 2 lines to load my data:
# SSE = load("SSE.jld")["SSE"]
# MSE = load("MSE.jld")["MSE"]
# If you run the 2 lines above, don't /5 later.
for σ in 𝛔
    σ_index = findfirst(𝛔 .== σ)
    for γ in 𝛄
        γ_index = findfirst(𝛄 .== γ)
        for i = 0:4
            # get 5-fold data:
            𝑿_test = permuted_matdata_train[(i*nrow_5fold + 1):((i+1) * nrow_5fold), 1:13]
            𝒚_test = permuted_matdata_train[(i*nrow_5fold + 1):((i+1) * nrow_5fold), 14]

            𝑿_train_1 = permuted_matdata_train[1:(i*nrow_5fold), 1:13]
            𝑿_train_2 = permuted_matdata_train[((1+i)*nrow_5fold+1):end, 1:13]
            𝑿_train = vcat(𝑿_train_1, 𝑿_train_2)

            𝒚_train_1 = permuted_matdata_train[1:(i*nrow_5fold), 14]
            𝒚_train_2 = permuted_matdata_train[((1+i)*nrow_5fold+1):end, 14]
            𝒚_train = vcat(𝒚_train_1, 𝒚_train_2)

            # compute 𝒚̂:
            𝑲 = zeros(𝑙, 𝑙)
            fill_𝑲(𝑲, 𝑿_train, σ)  # compute kernel matrix, # 𝐾ᵢⱼ = 𝐾(𝒙ᵢ, 𝒙ⱼ)

            𝜶 = get_𝜶(𝑲, γ, 𝑰, 𝒚_train)  # compute 𝛂*, 𝛂* = (𝑲 + γ𝑰)⁻¹⋅𝒚

            𝑿_test_vectorised = [𝑿_test[i,:] for i = 1:size(𝑿_test, 1)]
            get_ŷ_(x) = get_ŷ(𝜶, x, σ, γ, 𝑿_train)
            ŷ = get_ŷ_.(𝑿_test_vectorised)  # compute ŷ = 𝜶ᵀ⋅(𝐾₁ₐ, 𝐾₂ₐ, ..., 𝐾ₗₐ)ᵀ

            # have a look:
            display(plot(ŷ))
            display(plot!(𝒚_test))

            # compute testing error, MSE:
            sse = sum((𝒚_test - ŷ).^2)  # SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)²
            mse = sse/nrow_5fold  # MSE = SSE/N
            SSE[σ_index, γ_index] += sse
            MSE[σ_index, γ_index] += mse
        end
    end
end
SSE /= 5
MSE /= 5

# To save it so that I don't need to run again:
save("./data/SSE.jld", "SSE", SSE)
save("./data/MSE.jld", "MSE", MSE)

min_mse, σ_γ_indices = findmin(MSE)
σ_index, γ_index = σ_γ_indices[1], σ_γ_indices[2]

# Plot heatmap, take log so that the colourscale makes more sense:
heatmap(log2.(𝛄), log2.(𝛔), xticks=-40:-26, xlabel="gamma",
                            yticks=7:.5:13, ylabel="sigma", log.(log.(SSE)))
savefig("./graph/4.1.pdf")

############################## !! Comparing Everything !!
# 1 for naive regression
# 2-14 for 1-attribute linear regression
# 15 for all-attribute linear regression
jl4_𝐄 = zeros(16, 2)
jl4_𝛔 = zeros(16, 2)

# You can skip to last part by loading these:
# jl4_𝐄 = load("jl4_E.jld")["E"]
# jl4_𝛔 = load("jl4_s.jld")["s"]
############################## ADDING: Kernel Ridge Regression
# Run on whole training and testing sets:
σ = 𝛔[σ_index]
γ = 𝛄[γ_index]

𝑿_test  = permuted_matdata_test[:, 1:13]
𝑿_train = permuted_matdata_train[:, 1:13]
𝒚_test  = permuted_matdata_test[:, 14]
𝒚_train = permuted_matdata_train[:, 14]

𝑙 = nrow_train
𝑰 = Matrix{Float64}(I, 𝑙, 𝑙)  # 𝑰 = Id(l x l)

𝑲 = zeros(𝑙, 𝑙)
fill_𝑲(𝑲, 𝑿_train, σ)  # compute kernel matrix, # 𝐾ᵢⱼ = 𝐾(𝒙ᵢ, 𝒙ⱼ)

𝜶 = get_𝜶(𝑲, γ, 𝑰, 𝒚_train)  # compute 𝛂*, 𝛂* = (𝑲 + γ𝑰)⁻¹⋅𝒚

𝑿_train_vectorised = [𝑿_train[i,:] for i = 1:size(𝑿_train, 1)]
get_ŷ_(x) = get_ŷ(𝜶, x, σ, γ, 𝑿_train)
ŷ_train = get_ŷ_.(𝑿_train_vectorised)  # compute ŷ = 𝜶ᵀ⋅(𝐾₁ₐ, 𝐾₂ₐ, ..., 𝐾ₗₐ)ᵀ

𝑿_test_vectorised = [𝑿_test[i,:] for i = 1:size(𝑿_test, 1)]
get_ŷ_(x) = get_ŷ(𝜶, x, σ, γ, 𝑿_train)
ŷ_test = get_ŷ_.(𝑿_test_vectorised)  # compute ŷ = 𝜶ᵀ⋅(𝐾₁ₐ, 𝐾₂ₐ, ..., 𝐾ₗₐ)ᵀ

sse_train = sum((𝒚_train - ŷ_train).^2)  # SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)²
mse_train = sse_train/nrow_train  # MSE = SSE/N

sse_test = sum((𝒚_test - ŷ_test).^2)  # SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)²
mse_test = sse_test/nrow_test  # MSE = SSE/N

sse20_test  = zeros(20)
sse20_train = zeros(20)
for i = 1:20
    permuted_matdata = matdata[randperm(nrow),:]
    permuted_matdata_train = permuted_matdata[1:nrow_train,:]
    permuted_matdata_test = permuted_matdata[nrow_train+1:end,:]

    σ = 𝛔[σ_index]
    γ = 𝛄[γ_index]

    𝑿_test  = permuted_matdata_test[:, 1:13]
    𝑿_train = permuted_matdata_train[:, 1:13]
    𝒚_test  = permuted_matdata_test[:, 14]
    𝒚_train = permuted_matdata_train[:, 14]

    𝑙 = nrow_train
    𝑰 = Matrix{Float64}(I, 𝑙, 𝑙)  # 𝑰 = Id(l x l)

    𝑲 = zeros(𝑙, 𝑙)
    fill_𝑲(𝑲, 𝑿_train, σ)  # compute kernel matrix, # 𝐾ᵢⱼ = 𝐾(𝒙ᵢ, 𝒙ⱼ)

    𝜶 = get_𝜶(𝑲, γ, 𝑰, 𝒚_train)  # compute 𝛂*, 𝛂* = (𝑲 + γ𝑰)⁻¹⋅𝒚

    𝑿_train_vectorised = [𝑿_train[i,:] for i = 1:size(𝑿_train, 1)]
    get_ŷ_(x) = get_ŷ(𝜶, x, σ, γ, 𝑿_train)
    ŷ_train = get_ŷ_.(𝑿_train_vectorised)  # compute ŷ = 𝜶ᵀ⋅(𝐾₁ₐ, 𝐾₂ₐ, ..., 𝐾ₗₐ)ᵀ

    𝑿_test_vectorised = [𝑿_test[i,:] for i = 1:size(𝑿_test, 1)]
    get_ŷ_(x) = get_ŷ(𝜶, x, σ, γ, 𝑿_train)
    ŷ_test = get_ŷ_.(𝑿_test_vectorised)  # compute ŷ = 𝜶ᵀ⋅(𝐾₁ₐ, 𝐾₂ₐ, ..., 𝐾ₗₐ)ᵀ

    sse20_train[i] = sum((𝒚_train - ŷ_train).^2)  # SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)²
    #mse_train = sse_train/nrow_train  # MSE = SSE/N

    sse20_test[i] = sum((𝒚_test - ŷ_test).^2)  # SSE = 𝚺ᵢ(𝑦ᵢ - ̂𝑦ᵢ)²
    #mse_test = sse_test/nrow_test  # MSE = SSE/N
end

mse20_train = sse20_train/nrow_train
mse20_test  = sse20_test/nrow_test
plot(mse20_train, lab="training error")
plot!(mse20_test, lab="testing error")
savefig("./graph/4.2.pdf")

# x₁, ..., xₙ are n independent obeservations from a population
# that has mean μ and variance σ, then the variance of Σxᵢ is nσ²
# and the variance of x̄ = Σxᵢ/n is σ²/n
jl4_𝐄[16, 1] = 𝐄mse_train = mean(mse20_train)
jl4_𝐄[16, 2] = 𝐄mse_test  = mean(mse20_test)

jl4_𝛔[16, 1] = se_train = std(mse20_train)
jl4_𝛔[16, 2] = se_test  = std(mse20_test)

############################## ADDING: Naive Regression
# 20 runs of randomly splitted datasets:
ones_test  = ones(nrow_test)
ones_train = ones(nrow_train)
plot(matdata[:, 14])
sum_all_20_te  = zeros(20)  # no longer sum
sum_all_20_tse = zeros(20)
for i = 1:20
    # Obtain randomly splitted 𝒚_test/training:
    permuted_matdata = matdata[randperm(nrow),:]
    𝒚_train = permuted_matdata[nrow_test+1:nrow, 14]
    𝒚_test  = permuted_matdata[1:nrow_test, 14]

    # Plot each run:
    trl(x_test) = trained_regression_line(x_test, ones_train, 𝒚_train, 1)
    plot!(trl, 1, nrow)
    #println("done")

    # Obtain MSEs:
    sum_all_20_te[i]  = training_error_k_dim_basis(ones_train, 𝒚_train, 1)
    sum_all_20_tse[i] = test_error_k_dim_basis(ones_test, 𝒚_test, ones_train, 𝒚_train, 1)
end

# Show the result of plot loop by current():
current()

jl4_𝐄[1, 1] = 𝐄mse_naive_train = mean(sum_all_20_te)
jl4_𝐄[1, 2] = 𝐄mse_naive_test  = mean(sum_all_20_tse)

jl4_𝛔[1, 1] = se_naive_train = std(sum_all_20_te)
jl4_𝛔[1, 2] = se_naive_test  = std(sum_all_20_tse)

############################## ADDING: Linear Regression (attribute i)
# Using SuperLearn without basis function (but with integrated [ ,1]):
Core.eval(SuperLearn, :(TRANS_BASIS = false))

# the 𝑖th element is for the 𝑖th 𝑥:
sum_all_20_te  = zeros(14, 20)  # This is called sum for my historical reason
sum_all_20_tse = zeros(14, 20)
for j = 1:20
    permuted_matdata = matdata[randperm(nrow),:]
    for i = 1:13
        𝒙_test  = permuted_matdata[1:nrow_test, i]
        𝒙_train = permuted_matdata[nrow_test+1:nrow, i]
        𝒚_test  = permuted_matdata[1:nrow_test, 14]
        𝒚_train = permuted_matdata[nrow_test+1:nrow, 14]
        scatter(𝒙_train, 𝒚_train)
        trl(x_test) = trained_regression_line(x_test, 𝒙_train, 𝒚_train, nothing)
        display(plot!(trl, minimum(𝒙_train), maximum(𝒙_train)))

        sum_all_20_te[i, j]  = training_error_k_dim_basis(𝒙_train, 𝒚_train, nothing)
        sum_all_20_tse[i, j] = test_error_k_dim_basis(𝒙_test, 𝒚_test, 𝒙_train, 𝒚_train, nothing)
    end
end

jl4_𝐄[2:14, 1] = 𝐄mse_xi_train = [mean(sum_all_20_te[i,:])  for i in 1:13]
jl4_𝐄[2:14, 2] = 𝐄mse_xi_test  = [mean(sum_all_20_tse[i,:]) for i in 1:13]

jl4_𝛔[2:14, 1] = se_xi_train = [std(sum_all_20_te[i,:])  for i in 1:13]
jl4_𝛔[2:14, 2] = se_xi_test  = [std(sum_all_20_tse[i,:]) for i in 1:13]

############################## ADDING: Linear Regression (all attribute)
plot(matdata[:, 14])
X_test  = matdata[1:nrow_test, 1:13]
X_train = matdata[nrow_test+1:nrow, 1:13]
𝒚_test  = matdata[1:nrow_test, 14]
𝒚_train = matdata[nrow_test+1:nrow, 14]
trl(X_test) = trained_regression_line_M(X_test, X_train, 𝒚_train, nothing)
scatter!(nrow_test+1:nrow, trl(X_train))  # current() outside loop, display() inside!!

sorted_matdata = sort_matrix_by_jth_col(matdata, 14)
for j = 1:20
    permuted_matdata = matdata[randperm(nrow),:]
    X_test  = permuted_matdata[1:nrow_test, 1:13]
    X_train = permuted_matdata[nrow_test+1:nrow, 1:13]
    𝒚_test  = permuted_matdata[1:nrow_test, 14]
    𝒚_train = permuted_matdata[nrow_test+1:nrow, 14]

    sum_all_20_te[14, j]  = training_error_k_dim_basis(X_train, 𝒚_train, nothing)
    sum_all_20_tse[14, j] = test_error_k_dim_basis(X_test, 𝒚_test, X_train, 𝒚_train, nothing)
end

jl4_𝐄[15, 1] = 𝐄mse_xall_train = mean(sum_all_20_te[14,:])
jl4_𝐄[15, 2] = 𝐄mse_xall_test  = mean(sum_all_20_tse[14,:])

jl4_𝛔[15, 1] = se_xall_train = std(sum_all_20_te[14,:])
jl4_𝛔[15, 2] = se_xall_test  = std(sum_all_20_tse[14,:])

jl4_𝐄
jl4_𝛔
save("./data/jl4_E.jld", "E", jl4_𝐄)
save("./data/jl4_s.jld", "s", jl4_𝛔)

# Print Table:
# jl4_𝐄 = load("jl4_E.jld")["E"]
# jl4_𝛔 = load("jl4_s.jld")["s"]
jl4_𝐄_𝛔 = hcat(jl4_𝐄[:,1], jl4_𝛔[:,1], jl4_𝐄[:,2], jl4_𝛔[:,2])
jl4_𝐄_𝛔_rd = round.(jl4_𝐄_𝛔, digits=4)
jl4_𝐄_𝛔_str = (x -> @sprintf("%.2f", x)).(jl4_𝐄_𝛔_rd)
jl4_𝐄_str = (x -> @sprintf("%.2f", x)).(jl4_𝐄)
jl4_𝛔_str = (x -> @sprintf("%.2f", x)).(jl4_𝛔)
jl4 = jl4_𝐄_str .* " ± " .* jl4_𝛔_str

compare_table = convert(DataFrame, jl4)
names!(compare_table, [:E, :sigma])
compare_table[:Regression] = vcat(["Naive"], ["x$i" for i in 1:13], ["x"], ["KRR"])
compare_table = compare_table[[:Regression, :E, :sigma]]

# Save table to file:
print(compare_table)
io = open("./data/final_compare.txt", "w")
print(io, compare_table)
close(io)

# 16×3 DataFrame
# │ Row │ Regression │ E            │ sigma         │
# │     │ String     │ String       │ String        │
# ├─────┼────────────┼──────────────┼───────────────┤
# │ 1   │ Naive      │ 82.03 ± 6.00 │ 89.53 ± 12.30 │
# │ 2   │ x1         │ 72.15 ± 4.32 │ 71.45 ± 8.54  │
# │ 3   │ x2         │ 74.20 ± 4.12 │ 72.22 ± 8.27  │
# │ 4   │ x3         │ 65.28 ± 3.77 │ 63.75 ± 7.48  │
# │ 5   │ x4         │ 82.31 ± 4.31 │ 81.25 ± 8.66  │
# │ 6   │ x5         │ 69.85 ± 3.73 │ 67.59 ± 7.25  │
# │ 7   │ x6         │ 43.50 ± 4.35 │ 44.30 ± 9.29  │
# │ 8   │ x7         │ 73.26 ± 3.83 │ 71.03 ± 7.53  │
# │ 9   │ x8         │ 80.05 ± 4.53 │ 77.61 ± 9.00  │
# │ 10  │ x9         │ 72.69 ± 4.61 │ 71.30 ± 9.20  │
# │ 11  │ x10        │ 66.64 ± 4.22 │ 64.70 ± 8.38  │
# │ 12  │ x11        │ 62.96 ± 4.27 │ 62.41 ± 8.58  │
# │ 13  │ x12        │ 75.28 ± 4.30 │ 74.74 ± 8.49  │
# │ 14  │ x13        │ 38.70 ± 2.03 │ 38.32 ± 4.03  │
# │ 15  │ x          │ 22.03 ± 1.57 │ 23.13 ± 3.53  │
# │ 16  │ KRR        │ 8.13 ± 0.76  │ 12.90 ± 2.16  │

max_stars_repo_name: qin-yu/julia-regression-boston-housing

max_issues_repo_name: qin-yu/ml-julia-boston-housing

max_forks_repo_name: qin-yu/ml-julia-boston-housing


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

mutable struct ExpressRouteCircuitConnectionPropertiesFormat <: SwaggerModel
    expressRouteCircuitPeering::Any # spec type: Union{ Nothing, SubResource } # spec name: expressRouteCircuitPeering
    peerExpressRouteCircuitPeering::Any # spec type: Union{ Nothing, SubResource } # spec name: peerExpressRouteCircuitPeering
    addressPrefix::Any # spec type: Union{ Nothing, String } # spec name: addressPrefix
    authorizationKey::Any # spec type: Union{ Nothing, String } # spec name: authorizationKey
    ipv6CircuitConnectionConfig::Any # spec type: Union{ Nothing, Ipv6CircuitConnectionConfig } # spec name: ipv6CircuitConnectionConfig
    circuitConnectionStatus::Any # spec type: Union{ Nothing, CircuitConnectionStatus } # spec name: circuitConnectionStatus
    provisioningState::Any # spec type: Union{ Nothing, ProvisioningState } # spec name: provisioningState

    function ExpressRouteCircuitConnectionPropertiesFormat(;expressRouteCircuitPeering=nothing, peerExpressRouteCircuitPeering=nothing, addressPrefix=nothing, authorizationKey=nothing, ipv6CircuitConnectionConfig=nothing, circuitConnectionStatus=nothing, provisioningState=nothing)
        o = new()
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("expressRouteCircuitPeering"), expressRouteCircuitPeering)
        setfield!(o, Symbol("expressRouteCircuitPeering"), expressRouteCircuitPeering)
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("peerExpressRouteCircuitPeering"), peerExpressRouteCircuitPeering)
        setfield!(o, Symbol("peerExpressRouteCircuitPeering"), peerExpressRouteCircuitPeering)
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("addressPrefix"), addressPrefix)
        setfield!(o, Symbol("addressPrefix"), addressPrefix)
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("authorizationKey"), authorizationKey)
        setfield!(o, Symbol("authorizationKey"), authorizationKey)
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("ipv6CircuitConnectionConfig"), ipv6CircuitConnectionConfig)
        setfield!(o, Symbol("ipv6CircuitConnectionConfig"), ipv6CircuitConnectionConfig)
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("circuitConnectionStatus"), circuitConnectionStatus)
        setfield!(o, Symbol("circuitConnectionStatus"), circuitConnectionStatus)
        validate_property(ExpressRouteCircuitConnectionPropertiesFormat, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        o
    end
end # type ExpressRouteCircuitConnectionPropertiesFormat

const _property_map_ExpressRouteCircuitConnectionPropertiesFormat = Dict{Symbol,Symbol}(Symbol("expressRouteCircuitPeering")=>Symbol("expressRouteCircuitPeering"), Symbol("peerExpressRouteCircuitPeering")=>Symbol("peerExpressRouteCircuitPeering"), Symbol("addressPrefix")=>Symbol("addressPrefix"), Symbol("authorizationKey")=>Symbol("authorizationKey"), Symbol("ipv6CircuitConnectionConfig")=>Symbol("ipv6CircuitConnectionConfig"), Symbol("circuitConnectionStatus")=>Symbol("circuitConnectionStatus"), Symbol("provisioningState")=>Symbol("provisioningState"))
const _property_types_ExpressRouteCircuitConnectionPropertiesFormat = Dict{Symbol,String}(Symbol("expressRouteCircuitPeering")=>"SubResource", Symbol("peerExpressRouteCircuitPeering")=>"SubResource", Symbol("addressPrefix")=>"String", Symbol("authorizationKey")=>"String", Symbol("ipv6CircuitConnectionConfig")=>"Ipv6CircuitConnectionConfig", Symbol("circuitConnectionStatus")=>"CircuitConnectionStatus", Symbol("provisioningState")=>"ProvisioningState")
Base.propertynames(::Type{ ExpressRouteCircuitConnectionPropertiesFormat }) = collect(keys(_property_map_ExpressRouteCircuitConnectionPropertiesFormat))
Swagger.property_type(::Type{ ExpressRouteCircuitConnectionPropertiesFormat }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ExpressRouteCircuitConnectionPropertiesFormat[name]))}
Swagger.field_name(::Type{ ExpressRouteCircuitConnectionPropertiesFormat }, property_name::Symbol) =  _property_map_ExpressRouteCircuitConnectionPropertiesFormat[property_name]

function check_required(o::ExpressRouteCircuitConnectionPropertiesFormat)
    true
end

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

ext: jl


text:
using Metalhead, Test
using Flux

@testset "MLPMixer" begin
	@testset for mode in [:small, :base, :large] # :huge]
		@testset for drop_path_rate in [0.0, 0.5]
			m = MLPMixer(mode; drop_path_rate)
			@test size(m(x_224)) == (1000, 1)
			@test gradtest(m, x_224)
			GC.safepoint()
			GC.gc()
		end
	end
end

@testset "ResMLP" begin
  @testset for mode in [:small, :base, :large] # :huge]
		@testset for drop_path_rate in [0.0, 0.5]
			m = ResMLP(mode; drop_path_rate)
			@test size(m(x_224)) == (1000, 1)
			@test gradtest(m, x_224)
			GC.safepoint()
			GC.gc()
		end
  end
end

@testset "gMLP" begin
    @testset for mode in [:small, :base, :large] # :huge]
		@testset for drop_path_rate in [0.0, 0.5]
			m = gMLP(mode; drop_path_rate)
			@test size(m(x_224)) == (1000, 1)
			@test gradtest(m, x_224)
			GC.safepoint()
			GC.gc()
		end
    end
end


text:
"""
    gaussian_kernel(x1::Array{Float64,1}, x2::Array{Float64,1}, cov_func_parameters::GaussianKernelHyperparameters)
Returns a covariance estimated with a gaussian kernel.
"""
function gaussian_kernel(x1::AbstractArray{T,1}, x2::AbstractArray{R,1}, cov_func_parameters::GaussianKernelHyperparameters) where T<:Real where R<:Real
    return cov_func_parameters.w_0 * exp(-0.5 * sum(((x1 .- x2) ./ cov_func_parameters.w_i) .^ 2))
end

"""
    marginal_gaussian_kernel(x1::Array{Float64}, x2::Array{Float64}, cov_func_parameters::GaussianKernelHyperparameters)
Returns a covariance estimated with a gaussian kernel. Also returns the marginal covariances (how does each covariance change by bumping each hyperparameter).
"""
function marginal_gaussian_kernel(x1::AbstractArray{T,1}, x2::AbstractArray{R,1}, cov_func_parameters::GaussianKernelHyperparameters) where T<:Real where R<:Real
    # This should take in two points in N-dimensional space and a hyperparameter vector and output the covariance between them.
    # In addition it should output a matrix with all of the marginal values.
    w_i = cov_func_parameters.w_i
    dim_of_cov_func_parameters = length(w_i) + 1
    properval = gaussian_kernel(x1, x2, cov_func_parameters)
    marginal_vals = Array{typeof(cov_func_parameters.w_0),1}(undef, dim_of_cov_func_parameters)
    marginal_vals[1]  = properval/cov_func_parameters.w_0
    for i in 2:(dim_of_cov_func_parameters)
        marginal_vals[i] =  ((x1[i-1] - x2[i-1])^2 / w_i[i-1]^3 ) * properval
    end
    return (cov = properval, marginal_covariances = marginal_vals)
end

"""
    K_matrix(X::Array{Float64,2}, cov_func::Function, cov_func_parameters::GaussianKernelHyperparameters, noise::Float64 = 0.0)
Returns a K_matrix together with marginal K matrices (marginal over each hyperparameter). THe cov_func should be a function
with a signature like that of gaussian_kernel.
"""
function K_matrix(X::AbstractArray{R,2}, cov_func::Function, cov_func_parameters::GaussianKernelHyperparameters{<:Real}, noise::Real = 0.0) where R<:Real
    NoObs = size(X)[1]
    KK = diagm(0 => ones(NoObs))
    for r in 1:NoObs
        for c in (r+1):NoObs
            KK[r,c] = cov_func(X[r,:], X[c,:], cov_func_parameters)
        end
    end
    noise_matrix = noise * diagm(0 => ones(NoObs))
    return Hermitian(KK + noise_matrix)
end

"""
    correlation_vector_of_a_point(x::AbstractArray{T,1},X::AbstractArray{R,2}, cov_func::Function, cov_func_parameters::GaussianKernelHyperparameters) where T<:Real where R<:Real
This calculates the correlations of a point with each point in the X array.
"""
function correlation_vector_of_a_point(x::AbstractArray{T,1},X::AbstractArray{R,2}, cov_func::Function, cov_func_parameters::GaussianKernelHyperparameters) where T<:Real where R<:Real
    NoObs = size(X)[1]
    r = Array{R,1}(undef,NoObs )
    for i in 1:NoObs
        r[i] = cov_func(X[i,:], x, cov_func_parameters)
    end
    return r
end

"""
    K_matrix_with_marginals(X::Array{Float64,2}, cov_func::Function, cov_func_parameters::GaussianKernelHyperparameters, noise::Float64 = 0.0)
Returns a K_matrix together with marginal K matrices (marginal over each hyperparameter)
"""
function K_matrix_with_marginals(X::AbstractArray{T,2}, cov_func::Function, cov_func_parameters::GaussianKernelHyperparameters, noise::Real = 0.0) where T<:Real
    NoObs = size(X)[1]
    Ndims =  size(X)[2]
    number_of_marginal_matrices = length(cov_func_parameters.w_i) + 1
    covar_matrix = diagm(0 => ones(NoObs))
    mats = Array{typeof(cov_func_parameters.w_0),3}(undef, NoObs,NoObs,number_of_marginal_matrices)
    for m in 1:number_of_marginal_matrices
        mats[:,:,m] =  convert(Matrix{Float64}, zeros(Float64,NoObs,NoObs))
    end
    for r in 1:NoObs
        for c in (r+1):NoObs
            results, marginals = cov_func(X[r,:], X[c,:], cov_func_parameters)
            covar_matrix[r,c] = results
            for m in 1:number_of_marginal_matrices
                mats[r,c,m] = marginals[m]
            end
        end
    end
    cov_mat = Hermitian(covar_matrix + noise * diagm(0 => ones(NoObs)) )
    for m in 1:number_of_marginal_matrices
        mats[:,:,m] =  Hermitian(mats[:,:,m])
    end
    return (k_mat = cov_mat,  marginal_K_matrices = mats)
end

"""
    log_likelihood( y::Array{Float64,1},  K::Hermitian{Float64,Array{Float64,2}}; invK::Hermitian{Float64,Array{Float64,2}} = inv(K), determinant = det(K))
The log likelihood of a kriging model with values y and covariances K. invK and the determinant can be fed in as well to prevent additional operations.
Note that the normalising constant is excluded from the log likelihood here because it is not relevent for optimising hyperparameters.
"""
function log_likelihood(y::AbstractArray{T,1},  K::Hermitian; invK::Hermitian = inv(K), determinant::Real = det(K)) where T<:Real where R<:Real
    return -0.5 * transpose(y) * invK * y - 0.5 * log(determinant)
end

"""
    marginal_likelihood_gaussian_derivatives(X::Array{Float64,2}, y::Array{Float64,1}, w_0::Float64, w_i::Array{Float64,1}, noise::Float64 = 0.0)
The marginal likelihoods (along each parameter) of a kriging model are returned.
In addition the K matrix and the inverse K matrix are returned (to allow programers to use them as generated here and no redo them).
"""
function marginal_likelihood_gaussian_derivatives(X::AbstractArray{T,2}, y::AbstractArray{R,1}, cov_func_parameters::GaussianKernelHyperparameters, noise::Real = 0.0) where T<:Real where R<:Real
     w_0 = cov_func_parameters.w_0
     w_i = cov_func_parameters.w_i
    # From eqn 5.9 of Rasmussen and Williams
    K, marginal_covariances = K_matrix_with_marginals(X, marginal_gaussian_kernel, GaussianKernelHyperparameters(w_0,w_i), noise)
    invK = inv(K)
    alpha = invK * y
    alpha_alphaT__m_Kinv = alpha * transpose(alpha) - invK
    n_hyperparameters = size(marginal_covariances)[3]
    marginal_likelihoods = Array{Float64,1}(undef, n_hyperparameters)
    for i in 1:n_hyperparameters
        marginal_likelihoods[i] = 0.5 * tr(alpha_alphaT__m_Kinv * marginal_covariances[:,:,i])
    end
    return (marginal_likelihoods = marginal_likelihoods, k_mat = K, inv_K_matrix = invK)
end


text:
#  C:\dev\GitHub\autodiff\LowerTriangular.jl
#  AUTOGENERATED FROM C:\dev\GitHub\autodiff\LowerTriangular.ipynb on 2015-08-03T13:16:36
#[heading]
# Lower triangular matrices

using Base.Test
include("MatVec.jl")
# Make matrix from diagonal and strict lower triangle,
# e.g. D = [d11 d22 d33 d44]
#      LT = [L21 L31 L32 L41 L42 L43]
# Outputting
#  [d11   0   0   0]
#  [L21 d12   0   0] # row r: Ls starting at sum_i=1^r
#  [L31 L32 d33   0]
#  [L41 L42 L43 d44]
function ltri_unpack(D::Vec, LT::Vec)
  d=length(D)
  make_row(r::Int, L) = hcat(reshape([ L[i] for i=1:r-1 ],1,r-1), D[r], zeros(1,d-r))
  row_start(r::Int) = div((r-1)*(r-2),2)
  inds(r) = row_start(r)+(1:r-1)
  vcat([ make_row(r, LT[inds(r)]) for r=1:d ]...)
end

ltri_unpack(D, L) = ltri_unpack([Float64(d) for d in D], [Float64(l) for l in L])

@test [11 0 0 0; 21 22 0 0 ; 31 32 33 0 ; 41 42 43 44] == ltri_unpack([11 22 33 44], [21 31 32 41 42 43])

LL = ltri_unpack([1.1 2.2 3.3 4.4], [21 31 32 41 42 43])
@printf("An example lower triangle made from diag and LT=\n%s\n", LL)

function ltri_pack{M<:AbstractMatrix}(L::M)
  d=size(L,1)

  make_row(r::Int, L) = hcat(reshape([ L[i] for i=1:r-1 ],1,r-1), D[r], zeros(1,d-r))
  row_start(r::Int) = (r-1)*(r-2)/2
  diag(L), hcat([L[r,1:r-1] for r=1:d ]...)
end

ltri_pack{T}(L::LowerTriangular{T, Matrix{T}}) = ltri_pack(full(L)) ## Until  packed storage is implemented for those

@test ltri_unpack(ltri_pack(LL)...) == LL

@printf("packed=%s\n", ltri_pack(LL))

max_stars_repo_name: dsyme/ADBench

max_issues_repo_name: dsyme/ADBench

max_forks_repo_name: dsyme/ADBench


text:
struct SLArray{S,T,N,L,Syms} <: StaticArray{S,T,N}
  __x::SArray{S,T,N,L}
  #SLArray{Syms}(__x::StaticArray{S,T,N}) where {S,N,Syms,T} = new{S,N,Syms,T}(__x)
  SLArray{S,T,N,Syms}(__x::SArray) where {S,T,N,Syms} = new{S,T,N,length(__x),Syms}(convert.(T,__x))
  SLArray{S,Syms}(__x::SArray{S,T,N,L}) where {S,T,N,L,Syms} = new{S,T,N,L,Syms}(__x)
  SLArray{S,T,Syms}(__x::SArray{S,T,N,L}) where {S,T,N,L,Syms} = new{S,T,N,L,Syms}(__x)
  function SLArray{S,Syms}(x::Tuple) where {S,Syms}
    __x = SArray{S}(x)
    SLArray{S,Syms}(__x)
  end
  function SLArray{S,T,Syms}(x::Tuple) where {S,T,Syms}
    __x = SArray{S,T}(x)
    SLArray{S,T,Syms}(__x)
  end
  function SLArray{S,T,N,L,Syms}(x::Tuple) where {S,T,N,L,Syms}
    __x = SArray{S,T,N,L}(x)
    new{S,T,N,L,Syms}(__x)
  end
end

#####################################
# NamedTuple compatibility
#####################################
## SLArray to named tuple
function Base.convert(::Type{NamedTuple}, x::SLArray{S,T,N,L,Syms}) where {S,T,N,L,Syms}
  tup = NTuple{length(Syms),T}(x)
  NamedTuple{Syms,typeof(tup)}(tup)
end

## Named tuple to SLArray
#=
  1. `SLArray{Tuple{2,2}}((a=1, b=2, c=3, d=4))` (need to specify size)
  2. `SLArray{Tuple{2,2}}(a=1, b=2, c=3, d=4)` : alternative form using kwargs
  3. `SLVector((a=1, b=2))` : infer size for vectors
  4. `SLVector(a=1, b=2)` : alternative form using kwargs
=#
function SLArray{Size}(tup::NamedTuple{Syms,Tup}) where {Size,Syms,Tup}
  __x = Tup(tup) # drop symbols
  SLArray{Size,Syms}(__x)
end
SLArray{Size}(;kwargs...) where {Size} = SLArray{Size}(kwargs.data)
SLVector(tup::NamedTuple) = SLArray{Tuple{length(tup)}}(tup)
SLVector(;kwargs...) = SLVector(kwargs.data)

## pairs iterator
Base.pairs(x::SLArray{S,T,N,L,Syms}) where {S,T,N,L,Syms} =
    # (label => getproperty(x, label) for label in Syms) # not type stable?
    (Syms[i] => x[i] for i in 1:length(Syms))

#####################################
# StaticArray Interface
#####################################
@inline Base.getindex(x::SLArray, i::Int) = getfield(x,:__x)[i]
@inline Base.Tuple(x::SLArray) = Tuple(x.__x)

function StaticArrays.similar_type(::Type{SLArray{S,T,N,L,Syms}}, ::Type{NewElType},
    ::Size{NewSize}) where {S,T,N,L,Syms,NewElType,NewSize}
  n = prod(NewSize)
  if n == L
    SLArray{Tuple{NewSize...},NewElType,length(NewSize),L,Syms}
  else
    SArray{Tuple{NewSize...},NewElType,length(NewSize),n}
  end
end

function Base.similar(::Type{SLArray{S,T,N,L,Syms}}, ::Type{NewElType}, ::Size{NewSize}) where {S,T,N,L,Syms,NewElType,NewSize}
  n = prod(NewSize)
  if n == L
    tmp = Array{NewElType}(undef, NewSize)
    LArray{NewElType,length(NewSize),typeof(tmp),Syms}(tmp)
  else
    MArray{Tuple{NewSize...},NewElType,length(NewSize),n}(undef)
  end
end

Base.propertynames(::SLArray{S,T,N,L,Syms}) where {S,T,N,L,Syms} = Syms
symnames(::Type{SLArray{S,T,N,L,Syms}}) where {S,T,N,L,Syms} = Syms
@inline function Base.getproperty(x::SLArray,s::Symbol)
  s == :__x ? getfield(x,:__x) : x[s]
end
@inline function Base.getindex(x::SLArray,s::Symbol)
  syms = symnames(typeof(x))
  if syms isa NamedTuple
    idxs = syms[s]
    return idxs isa Tuple ? @views(x.__x[idxs...]) : @views(x.__x[idxs])
  else
    return getindex(x,Val(s))
  end
end
@inline @generated function Base.getindex(x::SLArray,::Val{s}) where s
    idx = findfirst(y->y==s,symnames(x))
    :(getfield(x,:__x)[$idx])
end

function Base.getindex(x::SLArray,inds::AbstractArray{I,1}) where I<:Integer
    getindex(x.__x,inds)
end
function Base.getindex(x::SLArray, inds::StaticVector{<:Any, Int})
    getindex(x.__x,inds)
end

# Note: This could in the future return an SLArray with the right names
function Base.getindex(x::SLArray,s::AbstractArray{Symbol,1})
    [getindex(x,si) for si in s]
end

"""
    @SLArray Size Names
    @SLArray Eltype Size Names

Creates an anonymous function that builds a labelled static vector with eltype
`ElType` with names determined from the `Names` with size `Size`. If no eltype
is given, then the eltype is determined from the arguments in the constructor.

For example:

```julia
ABC = @SLArray (2,2) (:a,:b,:c,:d)
x = ABC(1.0,2.5,3.0,5.0)
x.a == 1.0
x.b == 2.5
x.c == x[3]
x.d == x[2,2]
EFG = @SLArray (2,2) (e=1:3, f=4, g=2:4)
y = EFG(1.0,2.5,3.0,5.0)
EFG = @SLArray (2,2) (e=(2, :), f=4, g=2:4)
```

"""
macro SLArray(dims,syms)
  dims isa Expr && (dims = dims.args)
  syms = esc(syms)
  quote
    SLArray{Tuple{$dims...},$syms}
  end
end

macro SLArray(T,dims,syms)
  dims isa Expr && (dims = dims.args)
  syms = esc(syms)
  quote
    SLArray{Tuple{$dims...},$T,$(length(dims)),$(prod(dims)),$syms}
  end
end

"""
    @SLVector Names
    @SLVector Eltype Names

Creates an anonymous function that builds a labelled static vector with eltype
`ElementType` with names determined from the `Names`. If no eltype is given,
then the eltype is determined from the values in the constructor.

For example:

```julia
ABC = @SLVector (:a,:b,:c)
x = ABC(1.0,2.5,3.0)
x.a == 1.0
x.b == 2.5
x.c == x[3]
```

"""
macro SLVector(syms)
  syms = esc(syms)
  quote
    n = $syms isa NamedTuple ? maximum(map(maximum, $syms)) : length($syms)
    SLArray{Tuple{n},$syms}
  end
end

macro SLVector(T,syms)
  T = esc(T)
  syms = esc(syms)
  quote
    n = $syms isa NamedTuple ? maximum(map(maximum, $syms)) : length($syms)
    SLArray{Tuple{n},$T,1,n,$syms}
  end
end

"""
    symbols(::SLArray{T,N,D,Syms})

Returns the labels of the `SLArray` .

For example:

    z = SLVector(a=1, b=2, c=3)
    symbols(z)  # Tuple{Symbol,Symbol,Symbol} == (:a, :b, :c)
"""
symbols(::SLArray{S,T,N,L,Syms}) where {S,T,N,L,Syms} = Syms isa NamedTuple ? keys(Syms) : Syms


text:
# ------------------------------------------------------------------------------
# p-multigrid multilevel example
# ------------------------------------------------------------------------------

using LinearAlgebra

# setup
mesh = Mesh2D(1.0, 1.0)
finep = 4
midp = 2
coarsep = 1
numbercomponents = 1
dimension = 2
ctombasis =
    TensorH1LagrangePProlongationBasis(coarsep + 1, midp + 1, numbercomponents, dimension)
mtofbasis =
    TensorH1LagrangePProlongationBasis(midp + 1, finep + 1, numbercomponents, dimension)

# diffusion operators
finediffusion = GalleryOperator("diffusion", finep + 1, finep + 1, mesh)
middiffusion = GalleryOperator("diffusion", midp + 1, finep + 1, mesh)
coarsediffusion = GalleryOperator("diffusion", coarsep + 1, finep + 1, mesh)

# Chebyshev smoothers
finechebyshev = Chebyshev(finediffusion)
midchebyshev = Chebyshev(middiffusion)

# p-multigrid preconditioner
midmultigrid = PMultigrid(middiffusion, coarsediffusion, midchebyshev, [ctombasis])
multigrid = PMultigrid(finediffusion, midmultigrid, finechebyshev, [mtofbasis])

# compute operator symbols
A = computesymbols(multigrid, [3], [1, 1], [π, π])
eigenvalues = real(eigvals(A))

# ------------------------------------------------------------------------------


text:
@testset "Test add/remove services" begin
    sys = System(100)
    devices = []
    for i in 1:2
        bus = Bus(nothing)
        bus.name = "bus" * string(i)
        bus.number = i
        add_component!(sys, bus)
        gen = ThermalStandard(nothing)
        gen.bus = bus
        gen.name = "gen" * string(i)
        add_component!(sys, gen)
        push!(devices, gen)
    end

    service = StaticReserve{ReserveDown}(nothing)
    add_service!(sys, service, devices)

    for device in devices
        services = get_services(device)
        @test length(services) == 1
        @test services[1] isa Service
        @test services[1] == service
    end

    services = collect(get_components(Service, sys))
    @test length(services) == 1
    @test services[1] == service

    remove_component!(sys, service)

    for device in devices
        @test length(get_services(device)) == 0
    end
end

@testset "Test add_component Service" begin
    sys = System(100)
    static_reserve = StaticReserve{ReserveDown}(nothing)
    add_component!(sys, static_reserve)
    services = get_components(StaticReserve{ReserveDown}, sys)
    @test length(services) == 1
    @test iterate(services)[1] == static_reserve
end

@testset "Test add_service errors" begin
    sys = System(100)
    bus = Bus(nothing)
    service = StaticReserve{ReserveDown}(nothing)
    # Bus is not a Device.
    @test_throws ArgumentError add_service!(sys, service, [bus])

    gen = ThermalStandard(nothing)
    # gen is not in sys.
    @test_throws ArgumentError add_service!(sys, service, [bus])
end

@testset "Test remove service from device" begin
    sys = System(100)
    bus = Bus(nothing)
    bus.name = "bus1"
    bus.number = 1
    add_component!(sys, bus)
    gen = ThermalStandard(nothing)
    gen.bus = bus
    gen.name = "gen"
    add_component!(sys, gen)

    service = StaticReserve{ReserveDown}(nothing)
    add_service!(sys, service, [gen])
    @test length(get_services(gen)) == 1

    remove_service!(gen, service)
    @test length(get_services(gen)) == 0
end

@testset "Test has service" begin
    sys = System(100)
    bus = Bus(nothing)
    bus.name = "bus1"
    bus.number = 1
    add_component!(sys, bus)
    gen = ThermalStandard(nothing)
    gen.bus = bus
    gen.name = "gen"
    add_component!(sys, gen)

    service = StaticReserve{ReserveDown}(nothing)
    add_service!(sys, service, [gen])
    @test has_service(gen, service)
    @test has_service(gen, typeof(service))

    remove_service!(gen, service)
    @test !has_service(gen, service)
    @test !has_service(gen, typeof(service))
end

@testset "Test remove device with service" begin
    sys = System(100)
    bus = Bus(nothing)
    bus.name = "bus1"
    bus.number = 1
    add_component!(sys, bus)
    gen = ThermalStandard(nothing)
    gen.bus = bus
    gen.name = "gen"
    add_component!(sys, gen)

    service = StaticReserve{ReserveDown}(nothing)
    add_service!(sys, service, [gen])
    @test length(get_services(gen)) == 1

    remove_component!(sys, gen)
    @test length(get_services(gen)) == 0
end

@testset "Test clear_services" begin
    gen = ThermalStandard(nothing)
    service = StaticReserve{ReserveDown}(nothing)
    add_service!(gen, service)
    @test length(get_services(gen)) == 1

    remove_service!(gen, service)
    @test length(get_services(gen)) == 0
end

@testset "Test add device with service" begin
    sys = System(100)
    bus = Bus(nothing)
    bus.name = "bus1"
    bus.number = 1
    add_component!(sys, bus)

    gen = ThermalStandard(nothing)
    gen.bus = bus
    gen.name = "gen"

    service = StaticReserve{ReserveDown}(nothing)
    add_service!(gen, service)
    @test length(get_services(gen)) == 1

    @test_throws ArgumentError add_component!(sys, gen)
end

@testset "Test get_contributing_devices" begin
    sys = System(100)
    devices = []
    services = []
    for i in 1:5
        bus = Bus(nothing)
        bus.name = "bus" * string(i)
        bus.number = i
        add_component!(sys, bus)
        gen = ThermalStandard(nothing)
        gen.bus = bus
        gen.name = "gen" * string(i)
        add_component!(sys, gen)
        push!(devices, gen)

        service = StaticReserve{ReserveUp}(nothing)
        service.name = "StaticReserve" * string(i)
        push!(services, service)
        add_component!(sys, service)
    end

    add_service!(devices[1], services[1])
    add_service!(devices[2], services[1])
    add_service!(devices[3], services[2])
    add_service!(devices[4], services[2])
    add_service!(devices[5], services[2])

    expected_contributing_devices1 = [devices[1], devices[2]]
    expected_contributing_devices2 = [devices[3], devices[4], devices[5]]

    contributing_devices1 = get_contributing_devices(sys, services[1])
    contributing_devices2 = get_contributing_devices(sys, services[2])

    # Order of contributing_devices isn't guaranteed, sort them to test.
    sort!(contributing_devices1, by = x -> get_name(x))
    sort!(contributing_devices2, by = x -> get_name(x))
    @test contributing_devices1 == expected_contributing_devices1
    @test contributing_devices2 == expected_contributing_devices2

    mapping = get_contributing_device_mapping(sys)
    @test length(mapping) == length(services)
    key1 = ServiceContributingDevicesKey((StaticReserve{ReserveUp}, get_name(services[1])))
    key2 = ServiceContributingDevicesKey((StaticReserve{ReserveUp}, get_name(services[2])))
    key3 = ServiceContributingDevicesKey((StaticReserve{ReserveUp}, get_name(services[3])))
    @test haskey(mapping, key1)
    @test haskey(mapping, key2)
    @test haskey(mapping, key3)
    @test length(mapping[key1].contributing_devices) == 2
    @test length(mapping[key2].contributing_devices) == 3
    @test length(mapping[key3].contributing_devices) == 0

    sort!(mapping[key1].contributing_devices, by = x -> get_name(x))
    sort!(mapping[key2].contributing_devices, by = x -> get_name(x))
    @test mapping[key1].contributing_devices == expected_contributing_devices1
    @test mapping[key2].contributing_devices == expected_contributing_devices2
end

@testset "Test get_component combinations" begin
    sys = System(100)
    reserves = (
        StaticReserve{ReserveUp}(nothing),
        StaticReserve{ReserveDown}(nothing),
        VariableReserve{ReserveUp}(nothing),
        VariableReserve{ReserveDown}(nothing),
    )

    for reserve in reserves
        add_component!(sys, reserve)
    end

    @test length(get_components(Service, sys)) == length(reserves)
    @test length(get_components(Reserve, sys)) == length(reserves)
    @test length(get_components(StaticReserve, sys)) == 2
    @test length(get_components(VariableReserve, sys)) == 2
    @test length(get_components(StaticReserve{ReserveUp}, sys)) == 1
    @test length(get_components(StaticReserve{ReserveDown}, sys)) == 1
    @test length(get_components(VariableReserve{ReserveUp}, sys)) == 1
    @test length(get_components(VariableReserve{ReserveDown}, sys)) == 1
end

@testset "Test struct type collections" begin
    concrete_types = IS.get_all_concrete_subtypes(Service)
    reserve_types = InteractiveUtils.subtypes(Reserve)
    reserve_parametric_types = InteractiveUtils.subtypes(ReserveDirection)

    actual_count = length(concrete_types)
    for reserve in reserve_types
        for parametric in reserve_parametric_types
            actual_count += 1
        end
    end

    @test 5 == actual_count
end


text:
function make_scenarios(system::PowerSystems.System, scenario_quantity::Int64, variance::Float64)
    scenario_forecasts = Vector{PowerSystems.ScenarioBased}()
    for data_first_step in PowerSystems.get_forecast_initial_times(system)
        for f in PowerSystems.get_component_forecasts(PowerSystems.RenewableDispatch, system, data_first_step)
            ts = TimeSeries.timestamp(PowerSystems.get_timeseries(f))
            scenarios = Matrix{Float64}(undef, length(ts), scenario_quantity)
            for (ix, v) in enumerate(TimeSeries.values(PowerSystems.get_timeseries(f)))
                x = Random.rand(Distributions.Truncated(Distributions.Normal(v, variance*v), -Inf, 1.0), scenario_quantity)
                scenarios[ix, :] = x'
            end
            scenario_forecast = PowerSystems.ScenarioBased(PowerSystems.get_component(f), "PowerScenarios", TimeSeries.TimeArray(ts, scenarios))
            push!(scenario_forecasts, scenario_forecast)
        end
    end
    return scenario_forecasts
end


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

mutable struct ConnectivityHop <: SwaggerModel
    type::Any # spec type: Union{ Nothing, String } # spec name: type
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    address::Any # spec type: Union{ Nothing, String } # spec name: address
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    nextHopIds::Any # spec type: Union{ Nothing, Vector{String} } # spec name: nextHopIds
    issues::Any # spec type: Union{ Nothing, Vector{ConnectivityIssue} } # spec name: issues

    function ConnectivityHop(;type=nothing, id=nothing, address=nothing, resourceId=nothing, nextHopIds=nothing, issues=nothing)
        o = new()
        validate_property(ConnectivityHop, Symbol("type"), type)
        setfield!(o, Symbol("type"), type)
        validate_property(ConnectivityHop, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(ConnectivityHop, Symbol("address"), address)
        setfield!(o, Symbol("address"), address)
        validate_property(ConnectivityHop, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(ConnectivityHop, Symbol("nextHopIds"), nextHopIds)
        setfield!(o, Symbol("nextHopIds"), nextHopIds)
        validate_property(ConnectivityHop, Symbol("issues"), issues)
        setfield!(o, Symbol("issues"), issues)
        o
    end
end # type ConnectivityHop

const _property_map_ConnectivityHop = Dict{Symbol,Symbol}(Symbol("type")=>Symbol("type"), Symbol("id")=>Symbol("id"), Symbol("address")=>Symbol("address"), Symbol("resourceId")=>Symbol("resourceId"), Symbol("nextHopIds")=>Symbol("nextHopIds"), Symbol("issues")=>Symbol("issues"))
const _property_types_ConnectivityHop = Dict{Symbol,String}(Symbol("type")=>"String", Symbol("id")=>"String", Symbol("address")=>"String", Symbol("resourceId")=>"String", Symbol("nextHopIds")=>"Vector{String}", Symbol("issues")=>"Vector{ConnectivityIssue}")
Base.propertynames(::Type{ ConnectivityHop }) = collect(keys(_property_map_ConnectivityHop))
Swagger.property_type(::Type{ ConnectivityHop }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_ConnectivityHop[name]))}
Swagger.field_name(::Type{ ConnectivityHop }, property_name::Symbol) =  _property_map_ConnectivityHop[property_name]

function check_required(o::ConnectivityHop)
    true
end

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


text:
module Equil
using LightXML, Printf
using LinearAlgebra, Roots, RowEchelon
include("IdealGas.jl")
include("Utils.jl")

function equil(input_file::AbstractString, lib_dir::AbstractString)
    #lib directory
    xmldoc = parse_file(input_file)
    xmlroot = root(xmldoc)

    gasphase = get_collection_from_xml(xmlroot,"gasphase")
    thermo_file = lib_dir*"therm.dat"
    thermo_obj = IdealGas.create_thermo(gasphase, thermo_file)
    all_species_thermo = thermo_obj.thermo_all
    #create elements from the species composition
    elements = Array{String,1}()
    for sp in all_species_thermo
        append!(elements,collect(keys(sp.composition)))
    end
    unique!(elements)
    mole_fracs  = get_molefraction_from_xml(xmlroot,thermo_obj.molwt,gasphase)
    local T = get_value_from_xml(xmlroot,"T")
    local p = get_value_from_xml(xmlroot,"p")

    #Total number of moles under the given conditions
    n_total = p/R/T
    #Number of moles of each component
    moles = n_total*mole_fracs
    #Order the species according to inlet moles
    order_species(gasphase,moles)
    #Save the ordered list for later use
    gasphase_in = copy(gasphase)
    moles_in = copy(moles)

    #reconstruct the thermo obj based on the ordered list
    thermo_obj = IdealGas.create_thermo(gasphase, thermo_file)
    all_species_thermo = thermo_obj.thermo_all

    #Create the stoichiometric coefficients for the indepdendent reactions
    stc = stoichiometric_coeffient_matrix(gasphase,moles,elements,all_species_thermo)
    stc = round.(stc,digits=2)
    #The above method removes the inerts from the list. Therefore store the inlet moles of inerts for final calculation
    inert_species = setdiff(gasphase_in,gasphase)
    inert_moles = Array{Float64,1}()
    for i in inert_species
        append!(inert_moles,moles_in[get_index(i,gasphase_in)])
    end
    #Create the thermo_obj for the new list of species (i.e without inerts)
    thermo_obj = IdealGas.create_thermo(gasphase, thermo_file)
    all_species_thermo = thermo_obj.thermo_all

    #Uncomment this for printing the independent reactions
    # print_reactions(stc,gasphase)

    #get the Gibb's free energy of all species to calculate Kp
    H_all = IdealGas.H_all(thermo_obj,T)
    S_all = IdealGas.S_all(thermo_obj,T)
    G_all = H_all - T*S_all

    # #Calculate the Del G for each reaction
    DelG = Array{Float64,1}()
    for i in 1:size(stc,1)
        push!(DelG,sum(stc[i,:] .* G_all))
    end
    molefracs = moles_in/sum(moles_in)
    println("\nInititial condition:\n")
    println("Species \t moles \t\t molefraction")
    for k in eachindex(gasphase_in)
        @printf("%10s \t %.4e \t %.4e \n", gasphase_in[k], moles_in[k], molefracs[k])
    end
    Kp = exp.(-DelG/R/T)
    # Minimize the Gibb's energy of the mixture
    minimize_G!(moles, stc, thermo_obj, Kp, T, p)
    #Account for the inert moles
    moles_all = copy(moles)
    append!(moles_all,inert_moles)
    molefracs = moles_all/sum(moles_all)
    gasphase_in = copy(gasphase)
    append!(gasphase_in,inert_species)
    println("\nEquilibrium composition @ T= $T K and p=$p Pa\n")
    println("Species \t moles \t\t molefraction")
    for k in eachindex(gasphase_in)
        @printf("%10s \t %.4e \t %.4e\n", gasphase_in[k], moles_all[k],molefracs[k])
    end

    return Symbol("Success")
end

"""
Bubble sort the species with the highest to lower order of moles
"""
function order_species(gasphase::Array{String}, moles::Array{Float64})
    for i in 1:length(gasphase)-1
        for j in i+1:length(gasphase)
            if moles[i] < moles[j]
                moles[i], moles[j] = moles[j], moles[i]
                gasphase[i], gasphase[j] = gasphase[j], gasphase[i]
            end
        end
    end
end

"""
This function creates the stoichiometric coefficients for
    linearly indepdendent reactions based on null space
"""
function stoichiometric_coeffient_matrix(gasphase, moles, elements, all_species_thermo)
    #constuct the B matrix
    B = zeros(length(gasphase),length(elements))
    for i in eachindex(all_species_thermo)
        for (k,v) in all_species_thermo[i].composition
            j = get_index(String(k),elements)
            B[i,j] = v
        end
    end
    B = B'
    # Eliminate the inert element from B
    nr = Array{Int64,1}()
    for i in 1:size(B,1)
        if count(e->e>0, B[i,:]) < 2
            append!(nr,i)
        end
    end
    ncount = 0
    for i in nr
        B = B[1:end .!=i-ncount,:]
        ncount += 1
    end
    # Eliminate the element constituting the inert from the elements list
    deleteat!(elements,nr)
    # Eliminate the empty columns (delete the inert)
    nc = Array{Int64,1}()
    for i in 1:size(B,2)
        if count(e->e>0,B[:,i]) < 1
            append!(nc,i)
        end
    end
    ncount = 0
    for i in nc
        B = B[:,1:end .!= i-ncount]
        ncount += 1
    end
    #The inert needs to be removed from the gasphase list and the moles as well
    deleteat!(gasphase,nc)
    deleteat!(moles,nc)

    if size(B,1) != size(B,2)
        n_null_space = length(gasphase)-rank(B)
        A = zeros(n_null_space, length(gasphase)-n_null_space)
        AI = Matrix(I,n_null_space, n_null_space)
        B = vcat(B,hcat(A,AI))
        C = inv(B)
        nsv = C[:,end-n_null_space+1:end]
        b = minimum(abs.(filter(e->e!=0,nsv)))
        nsv /= b
        return nsv'
    else
        B = rref(A)
        #check for rows containing only zero elements and augment with identity matrix
        #to the right with size equal to the number of rows containing only zero as elements
        n_null_space = 0
        for i in 1:size(B,1)
            if count(e->e!=0, B(i,:)) == 0
                B[i,i] = 1
                n_null_space += 1
            end
        end
        C = inv(B)
        nsv = C[:,end-n_null_space+1:end]
        b = minimum(abs.(filter(e->e!=0,nsv)))
        nsv /= b
        return nsv'
    end

end

function print_reactions(stc::Matrix{Float64}, gasphase::Array{String})
    println("\nLinearly indepdendent reactions:\n")
    #Negative coefficients are assumed to be reactants
    for i in 1:size(stc,1)
        rcts = ""
        prdts = ""
        for k in 1:length(gasphase)
            if stc[i,k] < 0 && stc[i,k] != 0
                rcts *= string(abs(stc[i,k]))*" "*gasphase[k]*" + "
            elseif stc[i,k] > 0 && stc[i,k] != 0
                prdts *= string(stc[i,k])*" "*gasphase[k]*" + "
            end
        end
        rcts = (strip(rcts))[1:end-1]
        prdts = (strip(prdts))[1:end-1]
        println(rcts * " > " * prdts)
    end

end

function minimize_G!(moles::Array{Float64},stc::Matrix, thermo_obj, Kp::Array{Float64} , T::Float64, p::Float64)
    minimized = false
    nr = size(stc,1)
    Gmin = true
    while Gmin
        mole_fracs = moles/sum(moles)
        G0 = IdealGas.Gmix(thermo_obj,T,p,mole_fracs)
        for i in 1:nr
            equilibrate_reaction!(moles,stc[i,:],Kp[i],p)
            #check if the Gibbs energy has decreased
        end
        mole_fracs = moles/sum(moles)
        G = IdealGas.Gmix(thermo_obj,T,p,mole_fracs)
        # println(G, "\t", G0, "\t", G-G0)
        Gmin = G-G0 < 1e-15 ? false : true
    end
    # mole_fracs = moles/sum(moles)
    # println(mole_fracs)
    # readline()
end

function equilibrate_reaction!(moles::Array{Float64}, rxn_stc::Array{Float64}, Kp::Float64, p::Float64 )
    # println("Equilibrating rxn: ", rxn )
    dn_min = -1e10
    dn_max = 1e10
    #The maximum possible extent of a reaction depends on the limiting concentration that is present
    for i in eachindex(rxn_stc)
        if rxn_stc[i] < 0
            dn_max = min(dn_max,-moles[i]/rxn_stc[i])
        end
        if rxn_stc[i] > 0
            dn_min = max(dn_min,-moles[i]/rxn_stc[i])
        end
    end

    local moles_updated = copy(moles)

    function root_func(x)
        n1 = (moles_updated .+ rxn_stc*x)
        mole_fracs = n1/sum(n1)
        n_pp = mole_fracs * (p/p_std)
        prdt_p = 1.0
        prdt_r = 1.0
        for k in eachindex(rxn_stc)
            if rxn_stc[k] < 0
                prdt_r *= n_pp[k]^abs(rxn_stc[k])
            elseif rxn_stc[k] > 0
                prdt_p *= n_pp[k]^rxn_stc[k]
            end
        end
        return (Kp*prdt_r - prdt_p)
    end
    dn = find_zero(root_func,(dn_max,dn_min),Bisection())
    #moles = moles_updated .+ rxn*dn
    moles .+= rxn_stc*dn

end

end


text:
module ExperienceAnalysis

using Dates

export exposure

abstract type ExposurePeriod end

struct Anniversary{T} <: ExposurePeriod
    pol_period::T
end

struct AnniversaryCalendar{T,U} <: ExposurePeriod
    pol_period::T
    cal_period::U
end

struct Calendar{U} <: ExposurePeriod
    cal_period::U
end

# make ExposurePeriod broadcastable so that you can broadcast
Base.Broadcast.broadcastable(ic::ExposurePeriod) = Ref(ic)

function next_exposure(from, to, period)

	return (from = from, to = min(from + period, to))

end

"""
    exposure(ExposurePeriod,from,to,continued_exposure=false)

Return an array of name tuples `(from=Date,to=Date)` of the exposure periods for the given `ExposurePeriod`s.

If `continued_exposure` is `true`, then the final `to` date will continue through the end of the final ExposurePeriod. This is useful if you want the decrement of interest is the cause of termination, because then you want a full exposure.

# Example

```julia
julia> using ExperienceAnalysis,Dates

julia> issue = Date(2016, 7, 4)
julia> termination = Date(2020, 1, 17)
julia> basis = ExperienceAnalysis.Anniversary(Year(1))

julia> exposure(basis, issue, termination)
4-element Array{NamedTuple{(:from, :to),Tuple{Date,Date}},1}:
 (from = Date("2016-07-04"), to = Date("2017-07-04"))
 (from = Date("2017-07-04"), to = Date("2018-07-04"))
 (from = Date("2018-07-04"), to = Date("2019-07-04"))
 (from = Date("2019-07-04"), to = Date("2020-01-17"))

"""
function exposure(p::Anniversary{T},from, to, continued_exposure=false) where {T}
    period = p.pol_period
	result = [next_exposure(from, to, period)]
	while result[end].to < to
		push!(
            result,
            next_exposure(result[end].to, to, period)
            )
    end

    if continued_exposure && (result[end].to == to)
        result[end] = (from=result[end].from,to=result[end].from + period)
    end

	return result
end

function exposure(p::AnniversaryCalendar{T,U},from, to, continued_exposure=false) where {T,U}

    period = min(p.cal_period, p.pol_period)

    next_pol_per = from + p.pol_period
    next_cal_per = ceil(from, p.cal_period)

    next_terminus = min(min(next_pol_per, next_cal_per), to)

    result = [next_exposure(from, next_terminus, period)]
    while result[end].to < to
        while result[end].to < next_terminus
            push!(
                result,
                next_exposure(result[end].to, next_terminus, period)
                )
        end
        if result[end].to >= next_pol_per
            next_pol_per = next_pol_per + p.pol_period
        end
        if result[end].to >= next_cal_per
            next_cal_per = ceil(next_cal_per + p.cal_period, p.cal_period)
        end

        next_terminus = min(min(next_pol_per, next_cal_per), to)
    end

    if continued_exposure && (result[end].to == to)
        result[end] = (from=result[end].from,to=min(next_pol_per, next_cal_per))
    end

	return result
end

function exposure(p::Calendar{U},from, to,continued_exposure=false) where {U}
    period = p.cal_period

    next_cal_per = ceil(from, p.cal_period)

    next_terminus = min(next_cal_per, to)

    result = [next_exposure(from, next_terminus, period)]
    while result[end].to < to
        while result[end].to < next_terminus
            push!(
                result,
                next_exposure(result[end].to, next_terminus, period)
                )
        end
        if result[end].to >= next_cal_per
            next_cal_per = ceil(next_cal_per + p.cal_period, p.cal_period)
        end

        next_terminus = min(next_cal_per, to)
    end

    if continued_exposure && (result[end].to == to)
        result[end] = (from=result[end].from,to=result[end].from + period)
    end

	return result
end

end

max_stars_repo_name: JuliaActuary/ExperienceAnalysis

max_issues_repo_name: JuliaActuary/ExperienceAnalysis

max_forks_repo_name: JuliaActuary/ExperienceAnalysis


text:
function linesearch(xold::Vector{Float64},fold::Float64,
    g::Vector{Float64},p::Vector{Float64},stpmax::Float64,
    state::String,Vs::Float64,vs::Float64,ts::Float64,zs::Float64,
    Vlv::Float64,zeta::Float64,Q::Float64,rho::Float64,beta::Float64,W2root::Float64,
    c0::Float64,Kvo::Float64,Kvc::Float64,Ks::Float64,E::Float64,V0::Float64,
    A0::Float64,Aann::Float64,leff::Float64,h::Float64,f::Function,J::Function,
    maxiter::Int16)
    alpha = 1e-4;
    tolx = 1e-7;
    check = false;
    # Vs = 100*cm3Tom3;
    # vs = system.branches.c0[1][end];
    # if system.heart.av.zeta[n+1] == 0
    #     zs = 1e-5;
    # else
    #     zs = system.heart.av.zeta[n+1];
    # end
    #
    # f = CVModule.fav;
    # J = CVModule.Jav;

    # scale if attempted step too large
    np = norm(p);
    # println(p)
    if np > stpmax
        p = p*stpmax/np;
    end
    # println(p)

    slope = dot(g,p);
    if slope >= 0
        print(g)
        print(p)
        error("Roundoff problem in line search.")
    end

    # compute minimum step fraction
    test = 0.;
    for i = 1:length(xold)
        temp = abs(p[i])/max(abs(xold[i]),1.);
        if temp > test
            test = temp;
        end
    end
    alamin = tolx/test;

    omega = 1;
    alam = omega; # start w/ desired fraction of Newton step
    N = 1;
    Plv = 0.;
    Pao = 0.;

    while N < maxiter
        x = xold+alam*p;
        Plv = E.*(x[2].*Vs .- V0);
        Pao = beta.*((2.*rho./beta).*((x[1].*vs .- W2root)./8 .+
            c0).^2 .- sqrt.(A0));
        if Plv > Pao
            state = "opening";
        else
            state = "closing";
            # println(state)
        end
        # println(x)
        # println(xold)
        # println(alam)
        # println(p)
        JJ = J(x,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,state);
        # println(JJ)
        D = diagm(maximum!(zeros(length(x)),abs.(JJ)).^-1);
        fvec = D*f(x,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
            A0,Aann,leff,h,state);
        # println(D)
        # println(fvec)
        fn = 0.5*dot(fvec,fvec);
        if alam < alamin
            x = xold;
            # println(alam)
            # println(alamin)
            # println(fn)
            # println(fold+alpha*alam*slope)
            println("Δx converged in proximal line search. Verify in Newton loop.")
            check = true;
            return fn,x,check
        elseif fn < fold+alpha*alam*slope
            # println(fn)
            # println(fold)
            # println(alpha*alam*slope)
            return fn,x,check
        else
            if alam == omega
                tmplam = -slope/(2*(fn-fold-slope));
            else
                rhs1 = fn-fold-alam*slope
                rhs2 = fn2-fold-alam2*slope
                a = (rhs1/alam^2-rhs2/alam2^2)/(alam-alam2);
                b = (-alam2*rhs1/alam^2+alam*rhs2/alam2^2)/(alam-alam2);
                if a == 0
                    tmplam = -slope/(2*b);
                else
                    disc = b^2-3*a*slope;
                    if disc < 0
                        tmplam = 0.5*alam;
                    elseif b <= 0
                        tmplam = (-b+sqrt(disc)/(3*a));
                    else
                        tmplam = -slope/(b+sqrt(disc));
                    end
                end
                if tmplam > 0.5*alam
                    tmplam = 0.5*alam;
                end
            end
        end
        alam2 = alam;
        fn2 = fn;
        alam = max(tmplam,0.1*alam);
        N+=1;
        if N == maxiter
            println(xn)
            println(f(xn,Vs,vs,ts,zs,Vlv,zeta,Q,rho,beta,W2root,c0,Kvo,Kvc,Ks,E,V0,
                A0,Aann,leff,h,state))
            error("Proximal line search iteration failed to converge.");
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
@osx_only ENV["LLDB_DEBUGSERVER_PATH"] =
  "/Applications/Xcode.app/Contents/SharedFrameworks/LLDB.framework/Resources/debugserver"
using Gallium, Cxx
dbg = debugger()
include(Pkg.dir("Gallium","src","lldbrepl.jl"))
if isdefined(Base,:active_repl)
  RunLLDBRepl(dbg)
  Gallium.createTargetREPL(dbg)
  Gallium.RunTargetREPL(dbg)
  # Step up Target C++ mode
  Cxx.addHeaderDir(Gallium.TargetClang,joinpath(JULIA_HOME,"../../src"); kind = C_System)
  Cxx.addHeaderDir(Gallium.TargetClang,joinpath(JULIA_HOME,"../../src/support"); kind = C_System)
  Cxx.addHeaderDir(Gallium.TargetClang,joinpath(JULIA_HOME,"../../usr/include"); kind = C_System)
  cxxparse(Gallium.TargetClang,"""#include "julia.h" """)
  # Do this after julia.h
  Cxx.register_booth(Gallium.TargetClang)
  cxxparse(Gallium.TargetClang,readall(joinpath(dirname(@__FILE__),"../src/boottarget.cpp")))
end
lldb_exec(dbg,"target create $(joinpath(JULIA_HOME,"julia"))")
lldb_exec(dbg,"process attach --pid $(ARGS[1])")
lldb_exec(dbg,"thread select 1")
lldb_exec(dbg,"settings append target.source-map . $(joinpath(JULIA_HOME,"../../base"))")
lldb_exec(dbg,"settings set target.process.optimization-warnings false")
lldb_exec(dbg,"b jl_throw")


text:
export default_PARSDMM_options

"""
Returns a set of default options for the PARSDMM solver
"""
function default_PARSDMM_options(options,TF)

  if     TF == Float64
    TI = Int64
  elseif TF == Float32
    TI = Int32
  end

  options.x_min_solver          = "CG_normal" #what algorithm to use for the x-minimization (CG appied to normal equations)
  options.maxit                 = TI(200)     #max number of PARSDMM iterations
  options.evol_rel_tol          = TF(1e-3)    #stop PARSDMM if ||x^k - X^{k-1}||_2 / || x^k || < options.evol_rel_tol AND options.feas_tol is reached
  options.feas_tol              = TF(5e-2)    #stop PARSDMM if the transform-domain relative feasibility error is < options.feas_tol AND options.evol_rel_tol is reached
  options.obj_tol               = TF(1e-3)    #optional stopping criterion for change in distance from point that we want to project
  options.rho_ini               = [TF(10.0)]  #initial values for the augmented-Lagrangian penalty parameters. One value in array or one value per constraint set in array
  options.rho_update_frequency  = TI(2)       #update augmented-Lagrangian penalty parameters and relaxation parameters every X number of PARSDMM iterations
  options.gamma_ini             = TF(1.0)     #initial value for all relaxation parameters (scalar)
  options.adjust_rho            = true        #adapt augmented-Lagrangian penalty parameters or not
  options.adjust_gamma          = true        #adapt relaxation parameters in PARSDMM
  options.adjust_feasibility_rho= true        #adapt augmented-Lagrangian penalty parameters based on constraint set feasibility errors (can be used in combination with options.adjust_rho)
  options.Blas_active           = true        #use direct BLAS calls, otherwise the code will use Julia loop-fusion where possible
  options.feasibility_only      = false       #drop distance term and solve a feasibility problem
  options.FL                    = TF          #type of Float: Float32 or Float64
  options.parallel              = false       #comput proximal mappings, multiplier updates, rho and gamma updates in parallel
  options.zero_ini_guess        = true        #zero initial guess for primal, auxilliary, and multipliers
  Minkowski                     = false       #the intersection of sets includes a Minkowski set
  return options
end

max_stars_repo_name: slimgroup/SetIntersectionProjection

max_issues_repo_name: slimgroup/SetIntersectionProjection

max_forks_repo_name: slimgroup/SetIntersectionProjection


text:
FullDim_rec() = error("Cannot infer dimension of polyhedron constructed from no element.")
function FullDim_rec(it::It, its::Union{Rep, It}...)
    if vectortype(eltype(it)) <: StaticArrays.SVector
        return FullDim(eltype(it))
    elseif isempty(it)
        return FullDim_rec(its...)
    else
        return FullDim(first(it))
    end
end

FullDim_rep(rep::Rep, other_reps::Union{Nothing, Rep}...) = FullDim(rep)
FullDim_rep(rep::Nothing, other_reps::Union{Nothing, Rep}...) = FullDim_rep(other_reps...)


text:
function extract_bindings(exprs::AbstractArray)
    bindings = Set{Union{Expr,Symbol}}()
    for ex in exprs
        if isa(ex, Expr)
            # Positional parameter
            if ex.head == :(::)
                push!(bindings, ex.args[2])

            # Optional parameter
            elseif ex.head == :kw && isa(ex.args[1], Expr)
                push!(bindings, ex.args[1].args[2])
                union!(bindings, extract_bindings(ex.args[2:2]))

            # Keyword parameters
            elseif ex.head == :parameters
                union!(bindings, extract_bindings(ex.args))

            # Varargs parameter
            elseif ex.head == :...
                union!(bindings, extract_bindings(ex.args))

            # Default values for optional or keyword parameters
            elseif ex.head == :call
                push!(bindings, ex.args[1])
                union!(bindings, extract_bindings(ex.args[2:end]))
            end
        end
    end
    return bindings
end

function binding_expr(t::Type)
    joinbinding(fullname(t.name.module)..., t.name.name)
end
function binding_expr(f::Function)
    m = Base.function_module(f, Tuple)
    joinbinding(fullname(m)..., Base.function_name(f))
end

function adjust_bindings(translations::Dict)
    new_trans = Dict()
    for (k, v) in translations
        new_trans[k] = binding_expr(v)
    end
    return new_trans
end

function absolute_binding!(expr::Expr, translations::Dict)
    # Positional parameter
    if expr.head == :(::)
        binding = expr.args[2]
        expr.args[2] = translations[binding]

    # Optional parameter
    elseif expr.head == :kw && isa(expr.args[1], Expr)
        binding = expr.args[1].args[2]
        expr.args[1].args[2] = translations[binding]
        absolute_binding!(expr.args[2:2], translations)

    # Keyword parameters
    elseif expr.head == :parameters
        absolute_binding!(expr.args, translations)

    # Varargs parameter
    elseif expr.head == :...
        absolute_binding!(expr.args, translations)

    # Default values for optional or keyword parameters
    elseif expr.head == :call
        binding = expr.args[1]
        expr.args[1] = translations[binding]
        absolute_binding!(expr.args[2:end], translations)
    end
    expr
end

function absolute_binding!(exprs::Array, translations::Dict)
    for expr in exprs
        if isa(expr, Expr)
           absolute_binding!(expr, translations)
        end
    end
    exprs
end

function joinbinding(symbols::Symbol...)
    result = symbols[1]
    for s in symbols[2:end]
        result = Expr(:., result, QuoteNode(s))
    end
    return result
end

function splitbinding(expr::Union{Expr,Symbol})
    parts = Symbol[]
    if isa(expr, Expr) && expr.head == :.
        append!(parts, splitbinding(expr.args[1]))
        tail = expr.args[2]
        push!(parts, isa(tail, QuoteNode) ? tail.value : tail)
    else
        push!(parts, expr)
    end
    return parts
end

function variable_name(expr::Expr)
    # Positional parameter, with type assertion
    if expr.head == :(::)
        name = expr.args[1]  # x::Integer

    # Optional parameter
    elseif expr.head == :kw
        if isa(expr.args[1], Symbol)
            name = expr.args[1]  # x=0
        else
            name = expr.args[1].args[1]  # x::Integer=0
        end

    else
        error("Unable to process expression")
    end
    return name
end
variable_name(sym::Symbol) = sym

function call_parameters(expr::Expr)
    keywords = Expr[]
    positional = Any[]
    expr.head == :call || error("expression is not a call")
    for expr in expr.args[2:end]
        if isa(expr, Expr)
            # Keyword parameters
            if expr.head == :parameters
                for ex in expr.args
                    name = variable_name(ex)
                    push!(keywords, Expr(:kw, name, name))
                end

            # Varargs parameter
            elseif expr.head == :...
                name = variable_name(expr.args[1])
                push!(positional, Expr(:..., name))
            else
                push!(positional, variable_name(expr))
            end
        else
            push!(positional, variable_name(expr))
        end
    end
    if !isempty(keywords)
        positional = vcat(Expr(:parameters, keywords...), positional)
    end
    return positional
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

bat_default(
    ::typeof(bat_integrated_autocorr_len),
    ::Val{:algorithm},
    ::Union{AbstractVector{<:Real},AbstractVectorOfSimilarVectors{<:Real}}
) = GeyerAutocorLen()

bat_default(
    ::typeof(bat_eff_sample_size),
    ::Val{:algorithm},
    ::AbstractVectorOfSimilarVectors{<:Real}
) = GeyerAutocorLen()

bat_default(
    ::typeof(bat_eff_sample_size),
    ::Val{:algorithm},
    ::DensitySampleVector

) = GeyerAutocorLen()


text:
push!(LOAD_PATH, "../src/")

using Documenter
using Arrhenius

makedocs(
    format = Documenter.HTML(),
    sitename = "Arrhenius.jl",
    modules  = [Arrhenius],
    pages = [
        "Home" => "index.md"
    ]
)

deploydocs(;
    repo = "github.com/SuXY15/Arrhenius.jl",
    target = "build",
    devbranch = "main",
    devurl = "dev",
)


text:
abstract type AbstractKernel{T} end

return_type(K::AbstractKernel{T},args...) where {T} = T
ambient_dimension(K::AbstractKernel) = ambient_dimension(K.op)

struct SingleLayerKernel{T,Op} <: AbstractKernel{T}
    op::Op
end
SingleLayerKernel{T}(op) where {T} = SingleLayerKernel{T,typeof(op)}(op)
SingleLayerKernel(op)              = SingleLayerKernel{default_kernel_eltype(op)}(op)
(op::SingleLayerKernel)(x,y)       = op(SVector(x),SVector(y))

struct DoubleLayerKernel{T,Op} <: AbstractKernel{T}
    op::Op
end
DoubleLayerKernel{T}(op) where {T} = DoubleLayerKernel{T,typeof(op)}(op)
DoubleLayerKernel(op)              = DoubleLayerKernel{default_kernel_eltype(op)}(op)
(op::DoubleLayerKernel)(x,y,ny)       = op(SVector(x),SVector(y),SVector(ny))

struct AdjointDoubleLayerKernel{T,Op} <: AbstractKernel{T}
    op::Op
end
AdjointDoubleLayerKernel{T}(op) where {T} = AdjointDoubleLayerKernel{T,typeof(op)}(op)
AdjointDoubleLayerKernel(op)              = AdjointDoubleLayerKernel{default_kernel_eltype(op)}(op)
(op::AdjointDoubleLayerKernel)(x,y,nx)       = op(SVector(x),SVector(y),SVector(nx))

struct HyperSingularKernel{T,Op} <: AbstractKernel{T}
    op::Op
end
HyperSingularKernel{T}(op) where {T} = HyperSingularKernel{T,typeof(op)}(op)
HyperSingularKernel(op)              = HyperSingularKernel{default_kernel_eltype(op)}(op)
(op::HyperSingularKernel)(x,y,nx,ny) = op(SVector(x),SVector(y),SVector(nx),SVector(ny))

# kernel trait do distinguish various integral operators and potentials. This
# helps in dispatch and in known the expected kernel signature (e.g. does it
# require the normal at point x? What about at the ponit y? etc)
struct SingleLayer end
struct DoubleLayer end
struct AdjointDoubleLayer end
struct HyperSingular end

kernel_type(::SingleLayerKernel)        = SingleLayer()
kernel_type(::DoubleLayerKernel)        = DoubleLayer()
kernel_type(::AdjointDoubleLayerKernel) = AdjointDoubleLayer()
kernel_type(::HyperSingularKernel)      = HyperSingular()

combined_field_coefficients(::SingleLayerKernel)        = (0,-1)
combined_field_coefficients(::DoubleLayerKernel)        = (1,0)
combined_field_coefficients(::AdjointDoubleLayerKernel) = (0,-1)
combined_field_coefficients(::HyperSingularKernel)      = (1,0)

# most of the types, you want to specifically define this method for a given
# kernel K. In the generic case, however, we rely on promote_op to infer the
# return type given the trait `kernel_type`
function return_type(K,X,Y)
    x_type = eltype(qnodes(X))
    y_type = eltype(qnodes(Y))
    if kernel_type(K) == SingleLayer()
        return Base.promote_op(K,x_type,y_type)
    elseif kernel_type(K) == DoubleLayer()
        νy_type = eltype(qnormals(Y))
        return Base.promote_op(K,x_type,y_type,νy_type)
    elseif kernel_type(K) == AdjointDoubleLayer()
        νx_type = eltype(qnormals(X))
        return Base.promote_op(K,x_type,y_type,νx_type)
    elseif kernel_type(K) == HyperSingular()
        νx_type = eltype(qnormals(X))
        νy_type = eltype(qnormals(Y))
        return Base.promote_op(K,x_type,y_type,νx_type,νy_type)
    else
        error("unknown `kernel_type` trait")
    end
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
function skew_matrix(a::Vector{T})::Matrix{T} where T<:Real
    m = [  0    -a[3]  a[2]
           a[3]  0    -a[1]
          -a[2]  a[1]  0    ]
    return m
end

function predict(nominal::NominalState{T}, error::ErrorState{T}, error_cov::AbstractMatrix{T}, measurement::Measurement{T}, noise::Noise{T}, dt::T) where T<:Real
    R = rotationmatrix(nominal.q)

    # Nominal state kinematics
    new_nominal = NominalState(
        nominal.p + nominal.v * dt + 0.5 * (R * (measurement.am - nominal.ab) + nominal.g) * dt^2,
        nominal.v + (R * (measurement.am - nominal.ab) + nominal.g) * dt^2,
        nominal.q * qrotation((measurement.ωm - nominal.ωb) * dt),
        nominal.ab,
        nominal.ωb,
        nominal.g
    )

    # Jacobian of the error state kinematics
    Fx = zeros(6 * 3, 6 * 3)
    Fx[1:3, 1:3]     = Matrix{T}(I, 3, 3)
    Fx[1:3, 4:6]     = dt * Matrix{T}(I, 3, 3)
    Fx[4:6, 4:6]     = Matrix{T}(I, 3, 3)
    Fx[4:6, 7:9]     = -R * skew_matrix(measurement.am - nominal.ab) * dt
    Fx[4:6, 10:12]   = -R * dt
    Fx[4:6, 16:18]   = dt * Matrix{T}(I, 3, 3)
    Fx[7:9, 7:9]     = rotationmatrix(qrotation(measurement.ωm - nominal.ωb) * dt)'
    Fx[7:9, 13:15]   = -dt * Matrix{T}(I, 3, 3)
    Fx[10:18, 10:18] = Matrix{T}(I, 9, 9)

    # Error-state kinematics
    # new_error = ErrorState(Fx * vec(error))
    # We don't need to compute this, as the mean of the error is
    # always zero, so this equation always returns zero.

    # Jacobian of the random impulses kinematics
    Fi = zeros(6 * 3, 4 * 3)
    Fi[4:15, :] = Matrix{T}(I, 12, 12)

    # Covariance matrices of the random impulses
    Qi = zeros(4 * 3, 4 * 3)
    Qi[1:3, 1:3]     = noise.σan * dt^2 * Matrix{T}(I, 3, 3) # velocity noise Vi (m^2/s^2)
    Qi[4:6, 4:6]     = noise.σωn * dt^2 * Matrix{T}(I, 3, 3) # orientation noise Θi (rad^2)
    Qi[7:9, 7:9]     = noise.σaw * dt * Matrix{T}(I, 3, 3) # accelero bias noise Ai (m^2/s^4)
    Qi[10:12, 10:12] = noise.σωw * dt * Matrix{T}(I, 3, 3) # gyro bias noise Ωi (rad^2/s^2)

    # Error-state covariance
    new_error_cov = Fx * error_cov * Fx' + Fi * Qi * Fi'

    return new_nominal, error, new_error_cov
end


text:
# Assumption: the product of all elements could be stored in a variable::T
function productarray(numbers::Array{T,1})::Array{T,1} where T <: Number
    all = reduce( *, numbers )
    map( x -> all÷x, numbers )
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
##### Moisture tendencies

#####
##### First order fluxes
#####

function flux(::Tracers{N}, ::Advect, atmos, args) where {N}
    @unpack state = args
    u = state.ρu / state.ρ
    return (state.tracers.ρχ .* u')'
end

#####
##### Second order fluxes
#####

function flux(::Tracers{N}, ::Diffusion, atmos, args) where {N}
    @unpack state, aux, diffusive = args
    @unpack D_t = args.precomputed.turbulence
    d_χ = (-D_t) * aux.tracers.δ_χ' .* diffusive.tracers.∇χ
    return d_χ * state.ρ
end


text:
"""
    compute_frames!(elements::Vector{UA}; parent=nothing)
        where UA<:Union{AbstractObject,AbstractAction}

Set elem.frames.frames to the computed frames for each elem in elements.
"""
function compute_frames!(
    elements::Vector{UA};
    parent = nothing,
    parent_counter = 0,
) where {UA<:Union{AbstractObject,AbstractAction}}
    available_subframes = typemin(Int):typemax(Int)
    if parent !== nothing
        last_frames = get_frames(parent)
        available_subframes = 1:length(get_frames(parent))
    else
        last_frames = nothing
    end
    is_first = true
    counter = 1
    for elem in elements
        # update (CURRENT/PREVIOUS)_(OBJECT/ACTION) and CURRENT_OBJECT_ACTION_TYPE
        if elem isa Object
            set_current_object(elem)
            empty!(CURRENT_ACTION)
        elseif elem isa Action
            set_current_action(elem)
        end
        set_current_action_type(elem)

        if last_frames === nothing && get_frames(elem) === nothing
            throw(ArgumentError("Frames need to be defined explicitly in the initial
                Object/Background or Action."))
        end
        if get_frames(elem) === nothing
            set_frames!(parent, elem, last_frames; is_first = is_first)
        end
        last_frames = get_frames(elem)
        if !(get_frames(elem) ⊆ available_subframes)
            @warn(
                "Action defined outside the frame range of the parent object.
          Action #$counter for Object #$parent_counter is defined for frames
          $(get_frames(elem)) but Object #$parent_counter exists only for $(available_subframes).
          (Info: Background is counted as Object #1)"
            )
        end
        is_first = false
        counter += 1
    end
end

"""
    get_current_setting()

Return the current setting of the current object
"""
function get_current_setting()
    object = CURRENT_OBJECT[1]
    return object.current_setting
end

"""
    get_interpolation(frames::UnitRange, frame)

Return a value between 0 and 1 which represents the relative `frame` inside `frames`.
"""
function get_interpolation(frames::UnitRange, frame)
    frame == last(frames) && return 1.0
    t = (frame - first(frames)) / (length(frames) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
end

"""
    get_interpolation(action::AbstractAction, frame)

Return the value of the `action.anim` Animation based on the relative frame given by
`get_interpolation(get_frames(action), frame)`
"""
function get_interpolation(action::AbstractAction, frame)
    t = get_interpolation(get_frames(action), frame)
    if !(action.anim.frames[end].t ≈ 1)
        @warn "Animations should be defined from 0.0 to 1.0"
    end
    return interpolation_to_transition_val(at(action.anim, t), action.transition)
end

"""
    interpolation_to_transition_val(interpolation_val, Transition)

Returns the transition value for the given `interpolation_val`.
If the interpolation value is already of the correct form it just gets returned.
Otherwise the Transition function like `get_position` is called and the interpolated value
is calculated.
"""
interpolation_to_transition_val(t, ::Nothing) = t
interpolation_to_transition_val(t::Point, trans::Translation) = t
interpolation_to_transition_val(t::Float64, trans::Rotation) = t
interpolation_to_transition_val(t::Scale, trans::Scaling) = t

function interpolation_to_transition_val(t, trans::Translation)
    # does interpolate between `to` and `from` and assumes we are at `from` already
    from = get_position(trans.from)
    to = get_position(trans.to)
    return t * (to - from)
end

function interpolation_to_transition_val(t, trans::Scaling)
    from = get_scale(trans.from)
    to = get_scale(trans.to)
    return from + t * (to - from)
end

function isapprox_discrete(val; atol = 1e-4)
    return isapprox(val, round(val); atol = atol)
end

function polywh(polygon::Vector{Vector{Point}})
    T = typeof(polygon[1][1].x)
    min_x = typemax(T)
    min_y = typemax(T)
    max_x = typemin(T)
    max_y = typemin(T)
    for poly in polygon
        for p in poly
            min_x = min(min_x, p.x)
            min_y = min(min_y, p.y)
            max_x = max(max_x, p.x)
            max_y = max(max_y, p.y)
        end
    end
    return max_x - min_x, max_y - min_y
end

function get_polypoint_at(points, t; pdist = polydistances(points))
    if t ≈ 0
        return points[1]
    end
    ind, surplus = nearestindex(pdist, t * pdist[end])

    nextind = mod1(ind + 1, length(points))
    overshootpoint = between(
        points[ind],
        points[nextind],
        surplus / distance(points[ind], points[nextind]),
    )
    return overshootpoint
end

"""
    set_previous_object(object::Object)

Set the `object` as `PREVIOUS_OBJECT`
"""
function set_previous_object(object::Object)
    if isempty(PREVIOUS_OBJECT)
        push!(PREVIOUS_OBJECT, object)
    else
        PREVIOUS_OBJECT[1] = object
    end
end

"""
    set_previous_action(action::Action)

Set the `action` as `PREVIOUS_ACTION`
"""
function set_previous_action(action::Action)
    if isempty(PREVIOUS_ACTION)
        push!(PREVIOUS_ACTION, action)
    else
        PREVIOUS_ACTION[1] = action
    end
end

"""
    set_current_object(object::Object)

Set the `object` as `CURRENT_OBJECT` and update `PREVIOUS_OBJECT`/`PREVIOUS_ACTION`
"""
function set_current_object(object::Object)
    update_previous_object_or_action()

    if isempty(CURRENT_OBJECT)
        push!(CURRENT_OBJECT, object)
    else
        CURRENT_OBJECT[1] = object
    end
end

"""
    set_current_action(action::Action)

Set the `action` as `CURRENT_ACTION` and update `PREVIOUS_OBJECT`/`PREVIOUS_ACTION`
"""
function set_current_action(action::Action)
    update_previous_object_or_action()

    if isempty(CURRENT_ACTION)
        push!(CURRENT_ACTION, action)
    else
        CURRENT_ACTION[1] = action
    end
end

"""
    update_previous_object_or_action()

Update the `PREVIOUS_OBJECT` or `PREVIOUS_ACTION` depending on whether the
last element was an object or an action. This is still saved in `CURRENT_OBJECT_ACTION_TYPE`.
"""
function update_previous_object_or_action()
    if !isempty(CURRENT_OBJECT_ACTION_TYPE)
        if CURRENT_OBJECT_ACTION_TYPE[1] == :Object
            !isempty(CURRENT_OBJECT) && set_previous_object(CURRENT_OBJECT[1])
        else
            !isempty(CURRENT_ACTION) && set_previous_action(CURRENT_ACTION[1])
        end
    end
end

"""
    set_current_action_type(t)

Set `CURRENT_OBJECT_ACTION_TYPE` to `:Object` or `:Action` depending
on the type of `t`.
"""
function set_current_action_type(t)
    type = :Object
    if t isa AbstractAction
        type = :Action
    end
    if isempty(CURRENT_OBJECT_ACTION_TYPE)
        push!(CURRENT_OBJECT_ACTION_TYPE, type)
    else
        CURRENT_OBJECT_ACTION_TYPE[1] = type
    end
end

"""
    empty_CURRENT_constants()

empty all `CURRENT_` constants like `CURRENT_OBJECT`
"""
function empty_CURRENT_constants()
    empty!(CURRENT_VIDEO)
    empty!(CURRENT_OBJECT)
    empty!(CURRENT_ACTION)
    empty!(PREVIOUS_OBJECT)
    empty!(PREVIOUS_ACTION)
    empty!(CURRENT_OBJECT_ACTION_TYPE)
end

"""
    interpolateable(x::AbstractVector)

Return the vector in a datatype that is interpolateable.
Currently only implemented is to change from `<:Integer` to `float`
"""
interpolateable(x::AbstractVector) = x
interpolateable(x::AbstractVector{<:Integer}) = float.(x)

"""
    _get_range(sizemargin, sizefrom)

For even `sizemargin` returns a range with values from `sizemargin ÷ 2 + 1` to `sizefrom - sizemargin ÷ 2`.
For odd `sizemargin` the left range value is increased by one.
"""
function _get_range(sizemargin, sizefrom)
    return if sizemargin == 0
        1:sizefrom
    elseif iseven(sizemargin)
        (sizemargin ÷ 2 + 1):(sizefrom - sizemargin ÷ 2)
    else
        (sizemargin ÷ 2 + 2):(sizefrom - sizemargin ÷ 2)
    end
end

"""
    crop(im, heightto, widthto)

Crops an imagematrix to size `(heightto, widthto)` around the center of the image.
"""
function crop(im, heightto, widthto)
    heightfrom, widthfrom = size(im)
    heightmargin, widthmargin = (heightfrom - heightto), (widthfrom - widthto)
    heightrange = _get_range(heightmargin, heightfrom)
    widthrange = _get_range(widthmargin, widthfrom)
    cropped_im = im[heightrange, widthrange]
    return cropped_im
end

"""
    _apply_and_reshape(func, im, template, args...)

Applies function `func` to imagematrix `im`. Afterwards reshapes the output
cropping it or padding it to the size of `template`.
All the `args` are passed to func after `im` as arguments.
"""
function _apply_and_reshape(func, im, template, args...)
    newim = func(im, args...)
    template_size = size(template)
    if any(size(newim) .> template_size)
        newim = crop(newim, template_size...)
    end
    if any(size(newim) .< template_size)
        newim, _ = paddedviews(mean(im), newim, template)
    end
    return newim
end

"""
    default_postprocess(frame_image, frame, frames)

Returns its first argument. Used as effectless default for keywords
argument `postprocess_frame` in render function.
"""
function default_postprocess(frame_image, frame, frames)
    frame_image
end

"""
    _rescale_if_needed(frame_image, rescale_factor)

Converts `frame_image` values to RGB and if `rescale_factor`
in `render` is different from 1.0 rescales the images matrix for
faster rendering.
"""
function _convert_and_rescale_if_needed(frame_image, rescale_factor)
    frame_image = convert.(RGB, frame_image)
    return if !isone(rescale_factor)
        new_size = trunc.(Int, size(frame_image) .* rescale_factor)
        imresize(frame_image, new_size)
    else
        frame_image
    end
end

"""
    _postprocess(args...;kwargs...)

This function is used to perform the postprocessing inside the render function.
It does mainly two things:
 - Adds a memoization with checks for the frames to render. The keyword argument
frames_memory, stores those imagematrices associated to frames that will reappear later
in the frames processing, so that when those will have to be rendered they can be taken
from here instead, whereas those that do not reappear are discarded. (the check is performed
after each frame on those predent in the frames_memory)

 - Applies `postprocess_frame` as provided to `render` on each frame.
"""
function _postprocess(
    video,
    objects,
    frame;
    layers,
    frames_memory,
    postprocess_frame,
    frame_template,
    filecounter,
    frames,
    rescale_factor,
)
    # Memoized frame_image
    frame_image = if haskey(frames_memory, frame)
        frames_memory[frame]
    else
        frame_image = if frame == first(frames)
            frame_template
        else
            frame_image = _convert_and_rescale_if_needed(
                get_javis_frame(video, objects, frame; layers = layers),
                rescale_factor,
            )
        end
    end

    # Memoization container updating with cleaning
    if !(frame in frames[(filecounter + 1):end]) && haskey(frames_memory, frame)
        delete!(frames_memory, frame)
    elseif (frame in frames[(filecounter + 1):end]) && !haskey(frames_memory, frame)
        frames_memory[frame] = frame_image
    end

    return _apply_and_reshape(
        postprocess_frame,
        frame_image,
        frame_template,
        filecounter,
        frames,
    )
end