text:
function azimuthal_to_planar_eigenvector(::Type{Acoustic{T,3}},eigvecs::Array{Complex{T}}) where T
    basis_order = size(eigvecs,1) - 1
    S = size(eigvecs,2)
    P = size(eigvecs,3)

    v = [
        (m == 0) ? eigvecs[l+1,s,p] : zero(Complex{T})
    for p = 1:P for s = 1:S for l = 0:basis_order for m = -l:l]

    return reshape(v,basisorder_to_basislength(Acoustic{T,3},basis_order),S,P)
end

# haven't implemented the azimuthal symmetry equivalent for 2D
azimuthal_to_planar_eigenvector(::Type{Acoustic{T,2}},eigvecs::Array{Complex{T}}) where T = eigvecs

function eigensystem(ω::T, medium::Acoustic{T,2}, species::Species{T,2}, ::AbstractPlanarSymmetry;
        basis_order::Int = 2,
        numberofparticles::Number = Inf,
        kws...) where {T<:AbstractFloat}

    scale_number_density = one(T)
    if numberofparticles > 1
        scale_number_density = one(T) - one(T) / numberofparticles
    end

    k = ω / medium.c
    sps = species
    S = length(sps)
    ho = basis_order

    t_matrices = get_t_matrices(medium, sps, ω, ho)

    len = (2ho+1) * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    as = [
        s1.exclusion_distance * outer_radius(s1) + s2.exclusion_distance * outer_radius(s2)
    for s1 in sps, s2 in sps]

    function M_component(keff,Ns,j,l,m,n)
        (n == m ? 1.0 : 0.0)*(j == l ? 1.0 : 0.0) - 2.0pi * scale_number_density * number_density(sps[l]) * t_matrices[l][m+ho+1,m+ho+1] * Ns[n-m + 2ho+1,j,l] / (k^2.0-keff^2.0)
    end

    # this matrix is needed to calculate the eigenvectors
    function MM(keff::Complex{T})
        Ns = [
            kernelN2D(m,k*as[s1,s2],keff*as[s1,s2])
        for m = (-2ho):2ho, s1 = 1:S, s2 = 1:S]

        ind2 = 1
        for l = 1:S for n = -ho:ho
            ind1 = 1
            for s = 1:S for m = -ho:ho
                MM_mat[ind1,ind2] = M_component(keff,Ns,s,l,m,n)
                ind1 += 1
        end end
            ind2 += 1
        end end
        return MM_mat
    end

    return MM
end

function eigensystem(ω::T, medium::Acoustic{T,3}, species::Species{T,3}, ::AbstractPlanarSymmetry;
        basis_order::Int = 2,
        direction_eff::Union{AbstractVector{T},AbstractVector{Complex{T}}} = [0.0,0.0,1.0],
        numberofparticles::Number = Inf,
        kws...) where {T<:AbstractFloat}

    scale_number_density = one(T)
    if numberofparticles > 1
        scale_number_density = one(T) - one(T) / numberofparticles
    end

    k = ω/medium.c
    S = length(species)
    ho = basis_order
    len = basisorder_to_basislength(Acoustic{T,3},ho) * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    rθφp = cartesian_to_radial_coordinates(direction_eff)
    # below we take the real part because our associated Legendre functions are only implemented for real arguments.
    if norm(real.(rθφp) - rθφp) / norm(rθφp) < (eps(T))^(1/3)
        rθφp = real.(rθφp)
    end

    Ys = spherical_harmonics(2ho, rθφp[2], rθφp[3]);
    lm_to_n = lm_to_spherical_harmonic_index

    t_matrices = get_t_matrices(medium, species, ω, ho)
    t_diags = diag.(t_matrices)
    baselen(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)

    as = [
        s1.exclusion_distance * outer_radius(s1) + s2.exclusion_distance * outer_radius(s2)
    for s1 in species, s2 in species]

    function M_component(keff::Complex{T},Ns::Array{Complex{T}},l::Int,m::Int,s1::Int,dl::Int,dm::Int,s2::Int)::Complex{T}
        (m == dm && l == dl && s1 == s2 ? one(Complex{T}) : zero(Complex{T})) +
        4pi * as[s1,s2] * scale_number_density * number_density(species[s2]) * t_diags[s1][baselen(l)] *
        sum(
            Complex{T}(im)^(-l1) * Ys[lm_to_n(l1,dm-m)] * Ns[l1+1,s1,s2] *
            gaunt_coefficient(dl,dm,l,m,l1,dm-m)
        for l1 in max(abs(dm-m),abs(dl-l)):(dl+l)) / (keff^2.0 - k^2.0)
    end

    function MM(keff::Complex{T})::Matrix{Complex{T}}
        Ns = [kernelN3D(l,k*as[s1,s2],keff*as[s1,s2]) for l = 0:2ho, s1 = 1:S, s2 = 1:S]
        ind2 = 1
        for s2 = 1:S, dl = 0:ho for dm = -dl:dl
            ind1 = 1
            for s1 = 1:S, l = 0:ho for m = -l:l
                MM_mat[ind1, ind2] = M_component(keff,Ns,l,m,s1,dl,dm,s2)
                ind1 += 1
            end end
            ind2 += 1
        end end
        return MM_mat
    end

    return MM
end

function eigensystem(ω::T, medium::Acoustic{T,3}, species::Species{T,3}, ::PlanarAzimuthalSymmetry{3};
        basis_order::Int = 2,
        numberofparticles::Number = Inf,
        kws...) where {T<:AbstractFloat}

    scale_number_density = one(T)
    if numberofparticles > 1
        scale_number_density = one(T) - one(T) / numberofparticles
    end

    k = ω/medium.c
    S = length(species)
    ho = basis_order
    len = (ho+1) * S
    MM_mat = Matrix{Complex{T}}(undef,len,len)

    t_matrices = get_t_matrices(medium, species, ω, ho)
    t_diags = diag.(t_matrices)
    baselen(order::Int) = basisorder_to_basislength(Acoustic{T,3},order)

    as = [
        s1.exclusion_distance * outer_radius(s1) + s2.exclusion_distance * outer_radius(s2)
    for s1 in species, s2 in species]

    function M_component(keff::Complex{T},Ns::Array{Complex{T}},l,s1,dl,s2)::Complex{T}
        (l == dl && s1 == s2 ? one(Complex{T}) : zero(Complex{T})) +
        4pi * as[s1,s2] * scale_number_density * number_density(species[s2]) * t_diags[s1][baselen(l)] *
        sum(
            Complex{T}(im)^(-l1) * sqrt((2*l1+1)/(4pi) ) * Ns[l1+1,s1,s2] *
            gaunt_coefficient(dl,0,l,0,l1,0)
        for l1 in abs(dl-l):(dl+l)) / (keff^2.0 - k^2.0)
    end

    function MM(keff::Complex{T})::Matrix{Complex{T}}
        Ns = [kernelN3D(l,k*as[s1,s2],keff*as[s1,s2]) for l = 0:2ho, s1 = 1:S, s2 = 1:S]
        ind2 = 1
        for s2 = 1:S, dl = 0:ho
            ind1 = 1
            for s1 = 1:S, l = 0:ho
                MM_mat[ind1, ind2] = M_component(keff,Ns,l,s1,dl,s2)
                ind1 += 1
            end
            ind2 += 1
        end
        return MM_mat
    end

    return MM
end


text:
import .CUDA

function Base.cconvert(::Type{MPIPtr}, buf::CUDA.CuArray{T}) where T
    Base.cconvert(CUDA.CuPtr{T}, buf) # returns DeviceBuffer
end

function Base.unsafe_convert(::Type{MPIPtr}, X::CUDA.CuArray{T}) where T
    reinterpret(MPIPtr, Base.unsafe_convert(CUDA.CuPtr{T}, X))
end

# only need to define this for strided arrays: all others can be handled by generic machinery
function Base.unsafe_convert(::Type{MPIPtr}, V::SubArray{T,N,P,I,true}) where {T,N,P<:CUDA.CuArray,I}
    X = parent(V)
    pX = Base.unsafe_convert(CUDA.CuPtr{T}, X)
    pV = pX + ((V.offset1 + V.stride1) - first(LinearIndices(X)))*sizeof(T)
    return reinterpret(MPIPtr, pV)
end

function Buffer(arr::CUDA.CuArray)
    Buffer(arr, Cint(length(arr)), Datatype(eltype(arr)))
end

ext: jl


text:
function y_prob(X1::Array{T,3}, X2::Matrix{T}, D::Matrix{T}, α::Vector{T}, Π::Matrix{T},
    ξ_all::Vector{T}) where T
    N, J, K1, K2, L = get_dim(X1, X2, D)
    @assert length(ξ_all) == J
    y_prob = zeros(T, N, J)
    for (i, d) in enumerate(eachrow(D))
        index = exp.(X2*(Π*d) .+ X1[i,:,:]*α .+ ξ_all) # J-by-1
        y_prob[i,:] = index ./ sum(index)
    end
    return y_prob
end

# TODO: need to add a y_prob for hetero preference y_prob

"""
Sample likelihood when consumers face homogeneous choice set, no unobservable preference heterogeneity
"""
function llk_sample_fullset(Y::Vector, X1::Array, X2::Matrix, D::Matrix,
    α::Vector{T}, Π::Matrix{T}, ξ::Vector{T}) where T

    N, J, K1, K2, L = get_dim(X1, X2, D)
    ξ_all = vcat(zero(eltype(ξ)), ξ)
    @assert length(ξ_all) == J

    llk = 0.0
    for (i, d) in enumerate(eachrow(D))
        index = exp.(X2*(Π*d) .+ X1[i,:,:]*α .+ ξ_all) # J-by-1
        y_prob = index ./ sum(index)
        llk += log(y_prob[Y[i]])
    end
    return llk
end

"""
Sample likelihood when consumers face different choice sets, no unobservable preference heterogeneity

C::Matrix, N-by-J dummy matrix, determines the choice set of each individual
    This can be actual choice set restrictions, or purely random choice sets determined to facilitate estimation
"""
function llk_sample_subset(Y::Vector, X1::Array, X2::Matrix, D::Matrix, C::Vector{Vector{Int64}},
    α::Vector{T}, Π::Matrix{T}, ξ::Vector{T}) where T

    N, J, K1, K2, L = get_dim(X1, X2, D)
    ξ_all = vcat(zero(eltype(ξ)), ξ)
    @assert length(ξ_all) == J

    llk = 0.0
    for (i, d) in enumerate(eachrow(D))
        choice_id = findfirst(isequal(Y[i]), C[i])
        index = @views exp.(X2[C[i],:]*(Π*d) .+ X1[i,C[i],:]*α .+ ξ_all[C[i]]) # J-by-1
        y_prob = index ./ sum(index)
        llk += log(y_prob[choice_id])
    end
    return llk
end

function llk_sample(Y::Vector, X1::Array, X2::Matrix, D::Matrix,
    α::Vector{T}, Π::Matrix{T}, ξ::Vector{T}; subset = nothing) where T

    if isnothing(subset)
        llk = llk_sample_fullset(Y, X1, X2, D, α, Π, ξ)
    else
        llk = llk_sample_subset(Y, X1, X2, D, subset,  α, Π, ξ)
    end
    return llk
end

"""
Sample likelihood when there are unobservable preference heterogeneity
"""
function llk_sample(Y::Vector, X1::Array, X2::Matrix, D::Matrix,
    α::Vector{T}, Π::Matrix{T}, ξ::Vector{T}, σ::Vector{T}, n_sim::Integer) where T

    N, J, K1, K2, L = get_dim(X1, X2, D)
    ξ_all = vcat(zero(eltype(ξ)), ξ)
    @assert length(ξ_all) == J

    llk = 0.0
    for (i, d) in enumerate(eachrow(D))
        y_prob = zeros(T, n_sim, J)
        for i_sim = 1:n_sim
            index = exp.(X2*(Π*d .+ randn(K).*σ) .+ X1[i,:,:]*α .+ ξ_all) # J-by-1
            y_prob[i_sim,:] = index ./ sum(index)
        end
        y_prob = mean(y_prob, dims=1)
        llk += log(y_prob[Y[i]])
    end
    return llk
end

"""
CCP approach for homogeneous preference

    Note that the CCP method will only work for large number of commodities.
    Also, we will need IV for ξ to control endogeneity
"""
function ccp(Y::Vector, X::Matrix)
    # Y: N-by-1
    # X: J-by-K
    J, K = size(X)
    P = proportionmap(Y)
    δ = [P[j]-P[1] for j in 2:J]
    DX = zeros(J-1,K)
    for j = 2:J-1
        DX[j,:] = X[j,:] .- X[1,:]
    end
    β_hat = (DX'*DX)\(DX'*δ)
    return β_hat
end

function opt_homo(Y, X1, X2, D, init_guess; subset = nothing)
    N, J, K1, K2, L = get_dim(X1, X2, D)
    if isnothing(subset)
        f = para -> -llk_sample_fullset(Y, X1, X2, D,
            para[1:K1], # α
            reshape(para[K1+1:K1+K2*L], K2, L), # Π
            para[K1+K2*L+1:end]) # ξ
            # take negative to minimize
    else
        f = para -> -llk_sample_subset(Y, X1, X2, D, subset,
            para[1:K1], # α
            reshape(para[K1+1:K1+K2*L], K2, L), # Π
            para[K1+K2*L+1:end]) # ξ
            # take negative to minimize
    end
    opt = Optim.optimize(f, init_guess, LBFGS(); autodiff=:forward)
    return opt
end

function opt_hetero(Y, X, init_guess, n_sim)
    J, K = size(X)
    N = length(Y)
    f = para -> -llk_sample(Y, X, para[1:K], para[K+1:K+J-1], para[K+J:end]; n_sim=n_sim) # take negative to minimize, ξ only J-1 dim
    # g = TwiceDifferentiable(f, init_guess; autodiff=:forward)
    opt = Optim.optimize(f, init_guess, LBFGS(); autodiff=:forward)
    return opt
end

function estDCModel!(model::DCModel; init_guess::Vector, n_sim=5, subset=nothing)
    if !model.hetero_preference
        model.optResults = opt_homo(model.Y, model.X1, model.X2, model.D, init_guess; subset=subset)
    else
        model.optResults = opt_hetero(model.Y, model.X, init_guess, n_sim)
    end
    return model
end


text:
module SonarEquations

using SpecialFunctions

"""
	SE = signal_excess_passive(SL, TL, NL, DI, DT)

Calculates the signal excess in a passive sonar scenario for source level `SL` (dB), transmission loss `TL` (dB), noise level `NL` (dB), directivity index `DI` (dB) and detection threshold `DT` (dB).
"""
function signal_excess_passive(SL, TL, NL, DI, DT)
	SE = SL - TL - NL + DI - DT
end

"""
	detection_index_gaussian(p_dtc, p_fal)

Calculates the detection index for a non-fluctuating signal in a Gaussian-distributed noise environment for a probability of detection `p_dtc` and probability of false alarm `p_fal`.
"""
function detection_index_gaussian(p_dtc, p_fal)
	if p_fal ≥ p_dtc
		d = 0
	else
		# d = 2(erfcinv(2p_dtc) - erfcinv(2p_fal))^2
		d = 2(erfinv(1 - 2p_dtc) - erfinv(1 - 2p_fal))^2
	end
	return d
end
# function detection_index_gaussian(p_dtc, p_fal)
# 	d = 2(erfcinv(2p_dtc) - erfcinv(2p_fal))^2
# 	# d = 2(erfinv(1 - 2p_dtc) - erfinv(1 - 2p_fal))^2
# end

# function detection_index_gaussian(SL, TL, NL, B, t)
# 	d = B*t*((SL - TL)/(B*NL))^2
# end
function detection_index_gaussian(SL, TL, NL, B, t)
	if TL ≥ SL
		d = 0
	else
		d = B*t*((SL - TL)/NL)^2
	end
	return d
end

"""

"""
function detection_index_exponential(p_dtc, p_fal)
	d = (log(p_fal)/log(p_dtc) - 1)^2
end

"""

"""
function detection_index_rayleigh(p_dtc, p_fal)
	d = π/(4 - π)*(sqrt(log(p_fal)/log(p_dtc)) - 1)^2
end

"""

"""
function detection_threshold(d, B, t)
	DT = 5log10(d*B/t)
end

"""
	probability_of_detection_gaussian(d, p_fal)

Calculates the probability of detection for a non-fluctuating signal in a Gaussian-distributed noise environment for detection index `d` and probability of false alarm `p_fal`.
"""
function probability_of_detection_gaussian(d, p_fal)
	p_dtc = erfc(erfcinv(2p_fal) - sqrt(d/2))/2
end

"""

"""
function probability_of_detection_exponential(d, p_fal)
	p_dtc = exp(log(p_fal)/(sqrt(d) + 1))
end

"""

"""
function probability_of_detection_rayleigh(d, p_fal)
	p_dtc = exp(log(p_fal)/(sqrt((4 - π)/π*d) + 1)^2)
end

end


text:
using Test
using AdvancedMH
using Random
using Distributions
using Random

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

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

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

    # Construct a DensityModel.
    model = DensityModel(density)

    # Set up our sampler with initial parameters.
    spl = MetropolisHastings([0.0, 0.0])

    # Sample from the posterior.
    chain = sample(model, spl, 100000; param_names=["μ", "σ"])

    # chn_mean ≈ dist_mean atol=atol_v
    @test mean(chain["μ"].value) ≈ 0.0 atol=0.1
    @test mean(chain["σ"].value) ≈ 1.0 atol=0.1
end


text:
using HTTP, JSON3, StructTypes

# The order in which you define the structs matters!

struct Geo
  lat::String
  lng::String
end

struct Address
    street::String
    suite::String
    city::String
    zipcode::String
    geo::Geo
end

struct Company
    name::String
    catchPhrase::String
    bs::String
end

struct User
  id::Int
  name::String
  username::String
  email::String
  address::Address
  phone::String
  website::String
  company::Company
end

StructTypes.StructType(::Type{Geo}) = StructTypes.Struct()
StructTypes.StructType(::Type{Address}) = StructTypes.Struct()
StructTypes.StructType(::Type{Company}) = StructTypes.Struct()
StructTypes.StructType(::Type{User}) = StructTypes.Struct()

function getusers(url)
    response = HTTP.request("GET", url)
    return JSON3.read(String(response.body), Array{User})
end

# println(getusers("https://jsonplaceholder.typicode.com/users"))

max_stars_repo_name: vladflore/consume-rest

max_issues_repo_name: vladflore/consume-rest

max_forks_repo_name: vladflore/consume-rest


text:
using JudiLing
using CSV
using Test
using DataFrames

@testset "display" begin
    latin =
        DataFrame(CSV.File(joinpath("data", "latin_mini.csv")))
    cue_obj = JudiLing.make_cue_matrix(
        latin,
        grams = 3,
        target_col = :Word,
        tokenized = false,
        keep_sep = false,
    )
    JudiLing.display_matrix(latin, :Word, cue_obj, cue_obj.C, :C)

    n_features = size(cue_obj.C, 2)
    S = JudiLing.make_S_matrix(
        latin,
        ["Lexeme"],
        ["Person", "Number", "Tense", "Voice", "Mood"],
        ncol = n_features,
    )
    JudiLing.display_matrix(latin, :Word, cue_obj, S, :S)

    G = JudiLing.make_transform_matrix(S, cue_obj.C)
    JudiLing.display_matrix(latin, :Word, cue_obj, G, :G)

    Chat = S * G
    JudiLing.display_matrix(latin, :Word, cue_obj, Chat, :Chat)

    F = JudiLing.make_transform_matrix(cue_obj.C, S)
    JudiLing.display_matrix(latin, :Word, cue_obj, F, :F)

    Shat = cue_obj.C * F
    JudiLing.display_matrix(latin, :Word, cue_obj, Shat, :Shat)

    A = cue_obj.A
    JudiLing.display_matrix(latin, :Word, cue_obj, A, :A)

    acc, R = JudiLing.eval_SC(Chat, cue_obj.C, R = true)
    JudiLing.display_matrix(latin, :Word, cue_obj, R, :R)
end


text:
# An implementation of DIOM for the solution of the square linear system Ax = b.
#
# This method is described in
#
# Y. Saad, Iterative methods for sparse linear systems.
# PWS Publishing Company, Boston, USA, 1996.
#
# Y. Saad, Practical use of some krylov subspace methods for solving indefinite and nonsymmetric linear systems.
# SIAM journal on scientific and statistical computing, 5(1), pp. 203--228, 1984.
#
# Alexis Montoison, <alexis.montoison@polymtl.ca>
# Montreal, September 2018.

export diom

"""Solve the consistent linear system Ax = b using direct incomplete orthogonalization method.

DIOM is similar to CG with partial reorthogonalization.

An advantage of DIOM is that nonsymmetric or symmetric indefinite or both nonsymmetric
and indefinite systems of linear equations can be handled by this single algorithm.

This implementation allows a right preconditioner M.
"""
function diom(A :: AbstractLinearOperator, b :: AbstractVector{T};
              M :: AbstractLinearOperator=opEye(size(A,1)),
              atol :: Float64=1.0e-8, rtol :: Float64=1.0e-6,
              itmax :: Int=0, memory :: Int=20, verbose :: Bool=false) where T <: Number

  m, n = size(A)
  m == n || error("System must be square")
  length(b) == m || error("Inconsistent problem size")
  verbose && @printf("DIOM: system of size %d\n", n)

  # Initial solution x₀ and residual r₀.
  x = zeros(T, n)
  x_old = copy(x)
  # Compute β.
  rNorm = @knrm2(n, b) # rNorm = ‖r₀‖₂
  rNorm == 0.0 && return x, SimpleStats(true, false, [rNorm], [], "x = 0 is a zero-residual solution")

  iter = 0
  itmax == 0 && (itmax = 2*n)

  rNorms = [rNorm;]
  ε = atol + rtol * rNorm
  verbose && @printf("%5d  %7.1e\n", iter, rNorm)

  # Set up workspace.
  mem = min(memory, itmax) # Memory.
  V = [zeros(n) for i = 1 : mem] # Preconditioned Krylov vectors, orthogonal basis for {r₀, AM⁻¹r₀, (AM⁻¹)²r₀, ..., (AM⁻¹)ᵐ⁻¹r₀}.
  P = [zeros(n) for i = 1 : mem] # Directions for x : Pₘ = Vₘ(Uₘ)⁻¹.
  H = zeros(mem+2)  # Last column of the band hessenberg matrix Hₘ = LₘUₘ.
  # Each column has at most mem + 1 nonzero elements. hᵢ.ₘ is stored as H[m-i+2].
  # m-i+2 represents the indice of the diagonal where hᵢ.ₘ is located.
  # In addition of that, the last column of Uₘ is stored in H.
  L = zeros(mem) # Last mem Pivots of Lₘ.
  p = BitArray(undef, mem) # Last mem permutations.

  # Initial ξ₁ and V₁.
  ξ = rNorm
  @. V[1] = b / rNorm

  # Stopping criterion.
  solved = rNorm ≤ ε
  tired = iter ≥ itmax
  status = "unknown"

  while !(solved || tired)

    # Update iteration index.
    iter = iter + 1

    # Set position in circulars stacks.
    pos = mod(iter-1, mem) + 1 # Position corresponding to pₘ and vₘ in circular stacks P and V.
    next_pos = mod(iter, mem) + 1 # Position corresponding to vₘ₊₁ in the circular stack V.

    # Incomplete Arnoldi procedure.
    z = M * V[pos] # Forms pₘ
    w = A * z # Forms vₘ₊₁
    for i = max(1, iter-mem+1) : iter
      ipos = mod(i-1, mem) + 1 # Position corresponding to vᵢ in the circular stack V.
      diag = iter - i + 2
      H[diag] = @kdot(n, w, V[ipos]) # hᵢ.ₘ = < A * M⁻¹vₘ , vᵢ >
      @kaxpy!(n, -H[diag], V[ipos], w) # w ← w - hᵢ.ₘ * vᵢ
    end
    # Compute hₘ₊₁.ₘ and vₘ₊₁.
    H[1] = @knrm2(n, w) # hₘ₊₁.ₘ = ‖vₘ₊₁‖₂
    if H[1] ≠ 0.0 # hₘ₊₁.ₘ = 0 ⇒ "lucky breakdown"
      @. V[next_pos] = w / H[1] # vₘ₊₁ = w / hₘ₊₁.ₘ
    end
    # It's possible that uₘ₋ₘₑₘ.ₘ ≠ 0 when m ≥ mem + 1
    if iter ≥ mem + 2
      H[mem+2] = 0.0 # hₘ₋ₘₑₘ.ₘ = 0
    end

    # Update the LU factorization with partial pivoting of H.
    # Compute the last column of Uₘ.
    if iter ≥ 2
      for i = max(2,iter-mem+1) : iter
        lpos = mod(i-1, mem) + 1 # Position corresponding to lᵢ.ᵢ₋₁ in the circular stack L.
        diag = iter - i + 2
        next_diag = diag + 1
        if p[lpos]
          # The rows i-1 and i are permuted.
          H[diag], H[next_diag] = H[next_diag], H[diag]
        end
        # uᵢ.ₘ ← hᵢ.ₘ - lᵢ.ᵢ₋₁ * uᵢ₋₁.ₘ
        H[diag] = H[diag] - L[lpos] * H[next_diag]
      end
      # Compute ξₘ the last component of zₘ = β(Lₘ)⁻¹e₁.
      if !p[pos] # p[pos] ⇒ ξₘ = ξₘ₋₁
        # ξₘ = -lₘ.ₘ₋₁ * ξₘ₋₁
        ξ = - L[pos] * ξ
      end
    end

    # Compute the direction pₘ, the last column of Pₘ = Vₘ(Uₘ)⁻¹.
    for i = max(1,iter-mem) : iter-1
      ipos = mod(i-1, mem) + 1 # Position corresponding to pᵢ in the circular stack P.
      diag = iter - i + 2
      if ipos == pos
        # pₐᵤₓ ← -hₘ₋ₘₑₘ.ₘ * pₘ₋ₘₑₘ
        @kscal!(n, -H[diag], P[pos])
      else
        # pₐᵤₓ ← pₐᵤₓ - hᵢ.ₘ * pᵢ
        @kaxpy!(n, -H[diag], P[ipos], P[pos])
      end
    end
    # pₐᵤₓ ← pₐᵤₓ + M⁻¹vₘ
    @kaxpy!(n, 1.0, z, P[pos])

    # Determine if interchange between hₘ₊₁.ₘ and uₘ.ₘ is needed and compute next pivot lₘ₊₁.ₘ.
    if abs(H[2]) < H[1]
      p[next_pos] = true
      # pₘ = pₐᵤₓ / hₘ₊₁.ₘ
      @. P[pos] = P[pos] / H[1]
      # lₘ₊₁.ₘ = uₘ.ₘ / hₘ₊₁.ₘ
      L[next_pos] = H[2] / H[1]
    else
      p[next_pos] = false
      # pₘ = pₐᵤₓ / uₘ.ₘ
      @. P[pos] = P[pos] / H[2]
      # lₘ₊₁.ₘ = hₘ₊₁.ₘ / uₘ.ₘ
      L[next_pos] = H[1] / H[2]
    end

    # Compute solution xₘ.
    if p[pos]
      # xₘ = xₘ₋ₙ + ξₘ₋ₙ * pₘ
      # x_old = xₘ₋ₙ, with m-n is the last iteration without permutation at the next step
      @. x = x_old + ξ * P[pos]
    else
      # xₘ = xₘ₋₁ + ξₘ * pₘ
      @kaxpy!(n, ξ, P[pos], x)
    end

    # Update x_old and residual norm.
    if !p[next_pos]
      copyto!(x_old, x)
      # ‖ b - Axₘ ‖₂ = hₘ₊₁.ₘ * |ξₘ / uₘ.ₘ| without pivoting
      rNorm = H[1] * abs(ξ / H[2])
    else
      # ‖ b - Axₘ ‖₂ = hₘ₊₁.ₘ * |ξₘ / hₘ₊₁.ₘ| = |ξₘ| with pivoting
      rNorm = abs(ξ)
    end
    push!(rNorms, rNorm)

    # Update stopping criterion.
    solved = rNorm ≤ ε
    tired = iter ≥ itmax
    verbose && @printf("%5d  %7.1e\n", iter, rNorm)
  end
  verbose && @printf("\n")

  status = tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol"
  stats = SimpleStats(solved, false, rNorms, T[], status)
  return (x, stats)
end


text:
"""
    AbstractCumulant
"""
abstract type AbstractCumulant end

function Base.get(cumulant::AbstractCumulant, state_tp1, action_tp1, preds_tp1)
    throw(DomainError("get(CumulantType, args...) not defined!"))
end

"""
    FeatureCumulant <: AbstractCumulant

Basic Cumulant which takes the value c_tp1 = s_tp1[idx] for 1<=idx<=length(s_tp1)
"""
struct FeatureCumulant <: AbstractCumulant
    idx::Int
end

# Base.get(cumulant::FeatureCumulant, state_tp1, action_tp1, preds_tp1) =
#     state_tp1[cumulant.idx]

Base.get(cumulant::FeatureCumulant; kwargs...) =
        kwargs[:state_tp1][cumulant.idx]

"""
    PredictionCumulant <: AbstractCumulant

Basic cumulant which takes the value c_tp1 = preds_tp1[idx]
"""
struct PredictionCumulant <: AbstractCumulant
    idx::Int
end

Base.get(cumulant::PredictionCumulant, state_tp1, action_tp1, preds_tp1) =
    preds_tp1[cumulant.idx]

"""
    ScaledCumulant{F<:Number, T<:AbstractCumulant} <: AbstractCumulant

A cumulant which scales another AbstractCumulant
"""
struct ScaledCumulant{F<:Number, T<:AbstractCumulant} <: AbstractCumulant
    scale::F
    cumulant::T
end

Base.get(cumulant::ScaledCumulant, state_tp1, action_tp1, preds_tp1) =
    cumulant.scale*get(cumulant.cumulant, state_tp1, action_tp1, preds_tp1)

"""
    FunctionalCumulant{F<:Function} <: AbstractCumulant

A cumulant that has a user defined function c_tp1 = f(state_tp1, action_tp1, preds_tp1)
"""
struct FunctionalCumulant{F<:Function} <: AbstractCumulant
    f::F
end
#
# Base.get(cumulant::FunctionalCumulant, state_tp1, action_tp1, preds_tp1) =
#     cumulant.f(state_tp1, action_tp1, preds_tp1)

Base.get(cumulant::FunctionalCumulant; kwargs...) =
    cumulant.f(;kwargs...)

struct RewardCumulant <: AbstractCumulant end

# Base.get(cumulant::RewardCumulant; kwargs...) = @show kwargs, kwargs[:reward]
Base.get(cumulant::RewardCumulant; kwargs...) = kwargs[:reward]

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
module UKrig

using SpecialFunctions: besselk, gamma
import DynamicPolynomials: @polyvar, monomials, MonomialVector
using FixedPolynomials
using LinearAlgebra
using Statistics
using Random
using StaticArrays

export Mnu, krig_Mnu, krig_Mnu_plus, loglike_Mnu, loglike_Mnu_plus
export Gnu, krig_Gnu, krig_Gnu_plus, loglike_Gnu, loglike_Gnu_plus

## ===================================================
# Matérn auto-cov

tν𝒦t(t,ν) = t^ν * besselk(ν, t)

function Mnu(t, ν)::Float64
	pt, pν, p0, p1 = promote(t, ν, Float64(0), Float64(1)) # note besselk always returns a Float64 apparently
	return (pt==p0) ? p1 : tν𝒦t(√(2pν)*pt,pν) * 2^(1-pν) / gamma(pν)
end

# the const on Gnu which matches Matern's principle irregular term
function scν_Mnu(ν)
	if floor(ν)==ν
		return (2 * (-1)^ν)   * (- (ν/2)^ν / gamma(ν) / gamma(ν+1))
	else
		return (π / sin(ν*π)) * (- (ν/2)^ν / gamma(ν) / gamma(ν+1))
	end
end

## ===================================================
# Generalized auto-cov

scν(ν) = (-1)^(floor(Int,ν)+1)

function Gnu(t::A, ν::B) where {A<:Real, B<:Real}
	C = promote_type(A,B)
	if t==A(0)
		return C(0)
	end
	if floor(ν)==ν
		return C(scν(ν) * t^(2ν) * log(t))
	else
		return C(scν(ν) * t^(2ν))
	end
end

## ===================================================
# Util

ℓ2_dist(x,y) = norm(x - y)

function distmat(xcols::NTuple{d,Vx}, ycols::NTuple{d,Vy}) where {d,T<:Real,Vx<:AbstractVector{T}, Vy<:AbstractVector{T}}
	nx = length(xcols[1])
	ny = length(ycols[1])
	D = fill(T(0), nx, ny)
	for (x,y) in zip(xcols,ycols)
		D .+= (x .- y').^2
	end
	return sqrt.(D)
end

## ===================================================
# monomial closures

function _generate_fpb(monos::MonomialVector{true}, b::Vector, ::Val{d}) where {d}
    poly_fpb = Polynomial{Float64}.(dot(monos,b))
    cnf_sv = config(poly_fpb, SVector{d}(Vector{Float64}(undef,d)))
    fpb(x::SVector{d,Q}) where {Q<:Real} = evaluate(poly_fpb, x, cnf_sv)
    fpb(x::Real...) = fpb(SVector(x))
    fpb(x::NTuple{d,Q}) where {Q<:Real} = fpb(SVector(x))
    return fpb
end

function _construct_monos_Fp(m::Int, ::Val{d}) where {d}
    @polyvar x[1:d]
    monos   = monomials(x,0:m)
    poly_fp = System(Polynomial{Float64}.(monos))
    cnf_sv  = JacobianConfig(poly_fp, SVector{d}(Vector{Float64}(undef,d)))
    Fp(x::SVector{d,Q}) where {Q<:Real} = evaluate(poly_fp, x, cnf_sv)
    Fp(x::Real...) = Fp(SVector(x))
    Fp(x::NTuple{d,Q}) where {Q<:Real} = Fp(SVector(x))
    return monos, Fp
end

function _construct_Fmat_monos(m::Int, xdata::Vararg{Vector{T},d}) where {d,T<:Real}
	n   = length(xdata[1])
	monos, Fp = _construct_monos_Fp(m, Val(d))
	FpVec  = Fp.(xdata...)
	Fmat   = reduce(hcat,permutedims(FpVec))
	return Fmat, monos
end

function construct_constrasts(xdata::NTuple{d,Vector{T}}, mxdata::NTuple{ma,Vector{T}}; musr::Int=0, ν=0.5) where {d,ma,T<:Real}
	m = max(musr, floor(Int, ν))

	## monomial covariates
	Fmat_pre, monos = _construct_Fmat_monos(m, xdata...)
	mp = length(monos)

	## additional user supplied covariates
	if ma > 0
		Fmat = vcat(Fmat_pre,transpose.(mxdata)...)
	else
		Fmat = Fmat_pre
	end

	return Fmat, monos
end

## ===================================================
# Kriging closures

# """

# ```
# krig = krig_Gnu(fdata::Vector{T}, xdata::Vector{T}...; musr=0, nu=0.5, σg=1.0, σe=1.0)
# ```

# Here is a short summary how `krig` can be called and vectorized where `n` denotes the number of spatial
# interpolation points and `d` denotes the spatial dimension.

# ```
# xvar = Tuple(rand(n) for i=1:d)
# xs_in_rows = hcat(xvar...)
# xzip   = zip(xvar...)
# xvtupe = xzip |> collect
# xsv    = [SVector(x) for x in zip(xvar...)]
# xad    = [reduce(vcat,x) for x in zip(xvar...)]

# xvar isa NTuple{d,Vector{Float64}}
# xs_in_rows isa Matrix{Float64}
# xzip isa Base.Iterators.Zip{NTuple{d,Vector{Float64}}}
# xvtupe isa Vector{NTuple{d,Float64}}
# xsv isa Vector{SVector{d,Float64}}
# xad isa Vector{Vector{Float64}}
# ```

# Evaluation at a single spatial point
# ```
# krig(xsv[1])
# krig(xad[1]...)
# ```

# Evaluation at a collection of spatial points
# ```
# krig.(xsv)
# krig.(SVector{d}.(xad))
# krig.(xvar...)
# krig.(xzip)
# krig.(xvtupe)
# krig.(eachcol(xs_in_rows)...)
# ```

# """

function krig_Gnu(fdata::Vector{T}, xdata::Vararg{Vector{T},d}; musr::Int=0, σg=1.0, σe=1.0, ν=0.5) where {d,T<:Real}
	krig_pre, bma = krig_Gnu_plus(fdata, xdata, (); musr=musr, σg=σg, σe=σe, ν=ν)
	return krig_pre
end

function krig_Mnu(fdata::Vector{T}, xdata::Vararg{Vector{T},d}; musr::Int=0, σs=1.0, σe=1.0, ρ=1.0, ν=0.5) where {d,T<:Real}
	krig_pre, bma = krig_Mnu_plus(fdata, xdata, (); musr=musr, σg=σg, σe=σe, ρ=ρ, ν=ν)
	return krig_pre
end

function krig_Gnu_plus(fdata::Vector{T}, xdata::NTuple{d,Vector{T}}, mxdata::NTuple{ma,Vector{T}}; musr::Int=0, σg=1.0, σe=1.0, ν=0.5) where {d,ma,T<:Real}

	## spatial monomial covariates
	Fmat, monos = construct_constrasts(xdata, mxdata; musr=musr, ν=ν)

	## generalized auto-cov (without the var mult σg²)
	dmat = distmat(xdata, xdata)
	Gmat  = Gnu.(dmat, ν)

	## Solve for Krigin coeffs
	n     = length(fdata)
	mtot  = size(Fmat,1)
	mmono = length(monos)
	Ξ   = [
		(σg^2).*Gmat .+ (σe^2).*I(n)  Fmat'
		Fmat               zeros(mtot, mtot)
	]
	cb = Ξ \ vcat(fdata, zeros(mtot))
	c  = cb[1:n]
	b  = cb[n+1:end]
	bmono = b[1:mmono]
	bma   = b[(mmono+ma):end]
	fpb = _generate_fpb(monos, bmono, Val(d))

	## Generate Krigin closure
	function krig_pre(x::SVector{d,Q}) where Q<:Real
		sqdist = fill(Q(0),n)
		for i = 1:d
			sqdist .+= (x[i] .- xdata[i]).^2
		end
		Kvec = (σg^2).*Gnu.(sqrt.(sqdist), ν)
		return dot(Kvec,c) + fpb(x)
	end
    krig_pre(x::Real...) = krig_pre(SVector(x))
    krig_pre(x::NTuple{d,Q}) where {Q<:Real} = krig_pre(SVector(x))

	return krig_pre, bma # bma is used
end

# TODO: needs testing
function krig_Mnu_plus(fdata::Vector{T}, xdata::NTuple{d,Vector{T}}, mxdata::NTuple{ma,Vector{T}}; musr::Int=0, σg=1.0, σe=1.0, ρ=1.0, ν=0.5) where {d,ma,T<:Real}

	## spatial monomial covariates
	Fmat, monos = construct_constrasts(xdata, mxdata; musr=musr, ν=ν)

	## generalized auto-cov (without the var mult σg²)
	dmat = distmat(xdata, xdata)
	Mmat  = Mnu.(dmat ./ ρ, ν)

	## Solve for Krigin coeffs
	n     = length(fdata)
	mtot  = size(Fmat,1)
	mmono = length(monos)
	Ξ   = [
		(σs^2).*Mmat .+ (σe^2).*I(n)  Fmat'
		Fmat               zeros(mtot, mtot)
	]
	cb = Ξ \ vcat(fdata, zeros(mtot))
	c  = cb[1:n]
	b  = cb[n+1:end]
	bmono = b[1:mmono]
	bma   = b[(mmono+ma):end]
	fpb = _generate_fpb(monos, bmono, Val(d))

	## Generate Krigin closure
	function krig_pre(x::SVector{d,Q}) where Q<:Real
		sqdist = fill(Q(0),n)
		for i = 1:d
			sqdist .+= (x[i] .- xdata[i]).^2
		end
		Kvec = (σs^2).*Mnu.(sqrt.(sqdist), ν)
		return dot(Kvec,c) + fpb(x)
	end
    krig_pre(x::Real...) = krig_pre(SVector(x))
    krig_pre(x::NTuple{d,Q}) where {Q<:Real} = krig_pre(SVector(x))

	return krig_pre, bma # bma is used
end

## ===================================================
# REML closures for each fixed ν

function loglike_Gnu(fdata::Vector{T}, xdata::Vararg{Vector{T},d}; musr::Int=0, ν=0.5) where {d,T<:Real}
	return loglike_Gnu_plus(fdata, xdata, (); musr=musr, ν=ν)
end

function loglike_Mnu(fdata::Vector{T}, xdata::Vararg{Vector{T},d}; musr::Int=0, ν=0.5) where {d,T<:Real}
	return loglike_Mnu_plus(fdata, xdata, (); musr=musr, ν=ν)
end

function loglike_Gnu_plus(fdata::Vector{T}, xdata::NTuple{d,Vector{T}}, mxdata::NTuple{ma,Vector{T}}; musr::Int=0, ν=0.5) where {d,ma,T<:Real}

	Fmat, monos = construct_constrasts(xdata, mxdata; musr=musr, ν=ν)

	## generalized auto-cov (without the var mult σg²)
	dmat = distmat(xdata, xdata)
	Gmat  = Gnu.(dmat, ν)
    n   = length(fdata)
    Mᵀ = nullspace(Fmat)
    M  = transpose(Mᵀ)
	Mfdata = M * fdata

	function loglike_inv(σe, σg)
		Σ = M * ((σg^2).*Gmat .+ (σe^2).*I(n)) * Mᵀ  |> Symmetric
		Σ⁻¹Mfdata = Σ \ Mfdata
		ll = - (Mfdata⋅Σ⁻¹Mfdata) / 2 - logdet(Σ) / 2
		return ll
	end

	function loglike_chol(σe, σg)
		Σ = M * ((σg^2).*Gmat .+ (σe^2).*I(n)) * Mᵀ  |> Symmetric
		cholΣ = cholesky(Σ, check=false)
		if !LinearAlgebra.issuccess(cholΣ)
			return -Inf
		end
		L = cholΣ.L
		L⁻¹Mfdata = L \ Mfdata
		ll = - (L⁻¹Mfdata⋅L⁻¹Mfdata) / 2 - sum(log,diag(L))
		return ll
	end

	function wgrad_loglike_inv(σe, σg)
		uΣe  = M * Mᵀ  |> Symmetric
		uΣg  = M * Gmat * Mᵀ  |> Symmetric
		Σ    = (σg^2) * uΣg + (σe^2) * uΣe
		∂σeΣ = (2σe) * uΣe
		∂σgΣ = (2σg) * uΣg
		L = cholesky(Σ).L
		L⁻¹Mfdata = L \ Mfdata
		Σ⁻¹Mfdata = transpose(L) \ L⁻¹Mfdata
		ll    = - dot(L⁻¹Mfdata, L⁻¹Mfdata) / 2 - sum(log,diag(L))
		∂σell =  dot(Σ⁻¹Mfdata, ∂σeΣ * Σ⁻¹Mfdata) / 2 - tr(Σ \ ∂σeΣ) / 2
		∂σgll =  dot(Σ⁻¹Mfdata, ∂σgΣ * Σ⁻¹Mfdata) / 2 - tr(Σ \ ∂σgΣ) / 2
		return ll, ∂σell, ∂σgll
	end

	return loglike_inv, loglike_chol, wgrad_loglike_inv
end

# TODO: needs testing
function loglike_Mnu_plus(fdata::Vector{T}, xdata::NTuple{d,Vector{T}}, mxdata::NTuple{ma,Vector{T}}; musr::Int=0, ν=0.5) where {d,ma,T<:Real}

	Fmat, monos = construct_constrasts(xdata, mxdata; musr=musr, ν=ν)

	dmat  = distmat(xdata, xdata)
    n  = length(fdata)
    Mᵀ = nullspace(Fmat)
    M  = transpose(Mᵀ)
	Mfdata = M * fdata

	function loglike_inv(σe, σs, ρ)
		Mmat  = Mnu.(dmat ./ ρ, ν)
		Σ = M * ((σs^2).*Mmat .+ (σe^2).*I(n)) * Mᵀ  |> Symmetric
		Σ⁻¹Mfdata = Σ \ Mfdata
		ll = - (Mfdata⋅Σ⁻¹Mfdata) / 2 - logdet(Σ) / 2
		return ll
	end

	function loglike_chol(σe, σs, ρ)
		Mmat  = Mnu.(dmat ./ ρ, ν)
		Σ = M * ((σs^2).*Mmat .+ (σe^2).*I(n)) * Mᵀ  |> Symmetric
		cholΣ = cholesky(Σ, check=false)
		if !LinearAlgebra.issuccess(cholΣ)
			return -Inf
		end
		L = cholΣ.L
		L⁻¹Mfdata = L \ Mfdata
		ll = - (L⁻¹Mfdata⋅L⁻¹Mfdata) / 2 - sum(log,diag(L))
		return ll
	end

	return loglike_inv, loglike_chol
end

end # module


text:
@testset "20.valid-parentheses.jl" begin
    @test is_valid_parentheses("()") == true
    @test is_valid_parentheses("()[]{}") == true
    @test is_valid_parentheses("(]") == false
    @test is_valid_parentheses("([)]") == false
    @test is_valid_parentheses("{[]}") == true
end


text:
using  Statistics,StatsBase,Distributions,DataFrames,CSV,Plots,Match,Zygote
pyplot()

# Helper functions:
include("routines.jl")

## TEST SET CONSTRUCTION AND ANALYSIS
df=CSV.read("JigsawPred.csv",DataFrame;delim=",",header=false)

dfm=Matrix(df)
Ntest=size(df,1)
c_pred=zeros(Int,Ntest,6) # multi-label problem with 6 labels
for n=1:Ntest
    for k = 1:6
        if dfm[n,k]>0.5
            c_pred[n,k]=1
        end
    end
end
p1_pred = 1.0*dfm # for this classifier, the predictions are essentially deterministic

df2=CSV.read("JigsawTrue.csv",DataFrame;delim=",",header=false)
dfm=Matrix(df2)
c_true=zeros(Int,Ntest,6) # multi-label problem with 6 labels
for n=1:Ntest
    for k = 1:6
        if dfm[n,k]>0.5
            c_true[n,k]=1
        end
    end
end

F1Test_micro_true,F1Test_macro_true,accuracy_individual,accuracy_joint = MultiMetrics(c_pred,c_true)

# Bernoulli Sampling approximation:
lambda=0.9
p1_true_guess = lambda*p1_pred .+ (1-lambda)*0.5 # guess of the true label on test set

# Define the parameters for the Optimal Sampling Distribution:
htilde,dummy = happroxFn_multi_micro(c_pred,p1_true_guess)

Mvals=LinRange(10,Ntest,10) # expected number of samples

percent=zeros(length(Mvals))

acc_mean=zeros(length(Mvals))
acc_lower=zeros(length(Mvals))
acc_upper=zeros(length(Mvals))

f1_micro_mean=zeros(length(Mvals))
f1_micro_lower=zeros(length(Mvals))
f1_micro_upper=zeros(length(Mvals))

f1_macro_mean=zeros(length(Mvals))
f1_macro_lower=zeros(length(Mvals))
f1_macro_upper=zeros(length(Mvals))

for experiment=1:length(Mvals)

    global M = Mvals[experiment] # average number of elements in each test set replicate
    global b = optimalb(htilde,M)
    global s = (rand(Ntest) .< b)
    global total_samples_used = sum(s)

    ToLabel=findall(s.==1)

    F1hat_micro_mean,F1hat_micro_var=FhatBS_micro(c_pred[ToLabel,:],c_true[ToLabel,:],b[ToLabel])
    F1hat_macro_mean,F1hat_macro_var=FhatBS_macro(c_pred[ToLabel,:],c_true[ToLabel,:],b[ToLabel])

    # convert this to confidence bounds by fitting this to say a beta distribution:

    global confidence_level=0.9

    println("-----------------------------------------")
    println("total samples used  = $total_samples_used")
    println("% of test data      = $(total_samples_used/Ntest)")
    percent[experiment]=total_samples_used/Ntest

    println("")

    println("true F1             = $(F1Test_micro_true)")
    local mn,lower,upper = confidenceBeta(F1hat_micro_mean,F1hat_micro_var,confidence_level)
    println("estimated F1        = $mn [$lower,$upper]")

    f1_micro_mean[experiment]=F1hat_micro_mean
    f1_micro_lower[experiment]=lower
    f1_micro_upper[experiment]=upper

    println("")

    println("true F1             = $(F1Test_macro_true)")
    local mn,lower,upper = confidenceBeta(F1hat_macro_mean,F1hat_macro_var,confidence_level)
    println("estimated F1        = $mn [$lower,$upper]")

    f1_macro_mean[experiment]=F1hat_macro_mean
    f1_macro_lower[experiment]=lower
    f1_macro_upper[experiment]=upper

end

mp=maximum(percent)

f1_true = F1Test_micro_true
ribbon_upper=f1_micro_upper - f1_micro_mean
ribbon_lower=f1_micro_mean - f1_micro_lower
plot(reuse=false,show=true)
plot_f1_micro=plot!(percent,f1_micro_mean,ribbon=(ribbon_lower,ribbon_upper),label="$(confidence_level) confidence interval",framestyle = :box,minorgrid=true)
plot!(plot_f1_micro)
plot!([0,mp],[f1_true,f1_true],label="true value",legend=:bottomright)#,ylims=[0.5*f1_true,1.5*f1_true])
title!("F1 Micro")
xlabel!("fraction of test data labelled")

f1_true = F1Test_macro_true
ribbon_upper=f1_macro_upper - f1_macro_mean
ribbon_lower=f1_macro_mean - f1_macro_lower
plot(reuse=false,show=true)
plot_f1_macro=plot!(percent,f1_macro_mean,ribbon=(ribbon_lower,ribbon_upper),label="$(confidence_level) confidence interval",framestyle = :box,minorgrid=true)
plot!(plot_f1_macro)
plot!([0,mp],[f1_true,f1_true],label="true value",legend=:bottomright)#,ylims=[0.5*f1_true,1.5*f1_true])
title!("F1 Macro")
xlabel!("fraction of test data labelled")

max_stars_repo_name: davidbarber/TestEstimationPublic

max_issues_repo_name: davidbarber/TestEstimationPublic

max_forks_repo_name: davidbarber/TestEstimationPublic


text:
function MarkPerimetrNew(r)
    side = West
    counter = 0
    B = Array{Int, 1}(undef, 1000)
    B, counter = movetothecorner(r, side, B)

    for side1 in 3:6
        markline(r,HorizonSide(side1 % 4))
    end

    returnfromthecorner(r, inverse(side), B, counter)
end

function inverse(side)
    return HorizonSide((Int(side) + 2) % 4)
end

function movements(r, side)
    while !isborder(r, side)
        move!(r, side)
    end
end

function movetothecorner(r, side, B)
    counter = 1
    while ((!isborder(r, side)) || (!isborder(r, HorizonSide((Int(side) + 1) % 4))))
        if !isborder(r, side)
            move!(r, side)
            B[counter] = 1
            counter += 1
        end
        if !isborder(r, HorizonSide((Int(side) + 1) % 4))
            move!(r, HorizonSide((Int(side) + 1) % 4))
            B[counter] = 2
            counter += 1
        end
    end
    return B, counter
end

function returnfromthecorner(r, side, B, counter)
    for i in 1:(counter - 1)
        if B[counter - i] == 1
            move!(r, side)
        end
        if B[counter - i] == 2
            move!(r, HorizonSide((Int(side) + 1) % 4))
        end
    end
end

function markline(r, side)
    putmarker!(r)
    while isborder(r, side) != true
        move!(r,side)
        putmarker!(r)
    end
end


text:
# This loop was causing all memory to be allocated. It stopped,

# but all I did was add a println and then remove it.
#for (typef, ptype) in ((:int32, :Int32), (:int64 ,:Int64), (:uint64 ,:Uint64))
for (typef, ptype) in ((:Int32, :Int32), (:Int64 ,:Int64), (:UInt64 ,:UInt64))
     @eval begin
         p = genprimes($typef(1000))
         @test length(p) == 168
         @test eltype(p) == $ptype
     end
end

a = 24238423
b = 73923403
c = 100000

@test countprimes(b) == ntcountprimes(b)
@test typeof(genprimes("10")[1]) == UInt64
@test countprimes(a,b) == ntcountprimes(a,b)
# Don't want to import or require primes, I suppose
#@test genprimes(c) == Base.primes(c)
@test typeof(nprimes("10")[1]) == UInt64
@test ntcountprimes("10^9") == 50847534

## FIXME: bitrot,  get 234057667276344607 == -234057667276344607
#@test countprimes(:(10^19)) == 234057667276344607
## FIXME: bitrot
#@test countprimes("10^20") == 2220819602560918840

##@test typeof(primelookup("2^63")) == (Int64,(Int128,Int128,Int128))
@test typeof(primelookup("2^63")) == Tuple{Int64,Tuple{Int128,Int128,Int128}}

@test apopcount(zeros(10)) == 0
@test apopcount([]) == 0
@test apopcount([typemax(UInt64)]) == 64
@test apopcount([convert(UInt64,true)]) == 1
@test apopcount([true]) == 0

# fix bug. asking for limits below first table entry.
@test countprimes(7) == 4
@test countprimes(2,7) == 3


text:
require("Options")

module Optim
    using OptionsMod
    using Distributions
    using Calculus

    function centroid(p::Matrix)
         reshape(mean(p, 2), size(p, 1))
    end

    import Base.assign,
           Base.dot,
           Base.length,
           Base.push!,
           Base.ref,
           Base.show,
           Base.getindex,
           Base.setindex!

    export curve_fit,
           estimate_errors,
           optimize,
           DifferentiableFunction,
           TwiceDifferentiableFunction

    # Types
    include("types.jl")

    # Grid Search
    include("grid_search.jl")

    # Line Search Methods
    include("backtracking_line_search.jl")

    # Gradient Descent
    include("gradient_descent.jl")

    # Conjugate gradient
    include("cgdescent.jl")

    # Newton and Quasi-Newton Methods
    include("newton.jl")
    include("bfgs.jl")
    include("l_bfgs.jl")

    # Constrained optimization
    include("fminbox.jl")
    include("nnls.jl")

    # trust region methods
    include("levenberg_marquardt.jl")

    # Heuristic Optimization Methods
    include("nelder_mead.jl")
    include("simulated_annealing.jl")

    # End-User Facing Wrapper Functions
    include("optimize.jl")
    include("curve_fit.jl")
end


text:
@testset "Domains" begin
  @testset "Curve" begin
    c = Curve([j for i in 1:3, j in 1:5])
    @test ndims(c) == 3
    @test npoints(c) == 5
    @test coordtype(c) == Int
    @test coordinates(c, 1) == [1, 1, 1]
    @test coordinates(c, 5) == [5, 5, 5]

    @test sprint(show, c) == "5 Curve{Int64,3}"
    @test sprint(show, MIME"text/plain"(), c) == "5 Curve{Int64,3}\n 1  2  3  4  5\n 1  2  3  4  5\n 1  2  3  4  5"

    if visualtests
      gr(size=(800,800))
      c1 = Curve([j for i in 1:1, j in 1:5])
      c2 = Curve([j for i in 1:2, j in 1:5])
      c3 = Curve([j for i in 1:3, j in 1:5])
      @plottest plot(c1,ms=4) joinpath(datadir,"Curve1D.png") !istravis
      @plottest plot(c2,ms=4) joinpath(datadir,"Curve2D.png") !istravis
      @plottest plot(c3,ms=4) joinpath(datadir,"Curve3D.png") !istravis
    end
  end

  @testset "PointSet" begin
    ps = PointSet([1. 0.; 0. 1.])
    @test npoints(ps) == 2
    @test coordinates(ps, 1) == [1., 0.]
    @test coordinates(ps, 2) == [0., 1.]

    @test sprint(show, ps) == "2 PointSet{Float64,2}"
    @test sprint(show, MIME"text/plain"(), ps) == "2 PointSet{Float64,2}\n 1.0  0.0\n 0.0  1.0"

    if visualtests
      gr(size=(800,800))
      Random.seed!(2019)
      @plottest plot(PointSet(rand(1,10))) joinpath(datadir,"pset1D.png") !istravis
      @plottest plot(PointSet(rand(2,10))) joinpath(datadir,"pset2D.png") !istravis
      @plottest plot(PointSet(rand(3,10))) joinpath(datadir,"pset3D.png") !istravis
      @plottest plot(PointSet(rand(1,10)),1:10) joinpath(datadir,"pset1D-data.png") !istravis
      @plottest plot(PointSet(rand(2,10)),1:10) joinpath(datadir,"pset2D-data.png") !istravis
      @plottest plot(PointSet(rand(3,10)),1:10) joinpath(datadir,"pset3D-data.png") !istravis
    end
  end

  @testset "RegularGrid" begin
    grid = RegularGrid{Float32}(200,100)
    @test ndims(grid) == 2
    @test coordtype(grid) == Float32
    @test size(grid) == (200,100)
    @test npoints(grid) == 200*100
    @test coordinates(grid, 1) == [0.,0.]
    @test origin(grid) == (0f0, 0f0)
    @test spacing(grid) == (1f0, 1f0)

    grid = RegularGrid((200,100,50), (0.,0.,0.), (1.,1.,1.))
    @test ndims(grid) == 3
    @test coordtype(grid) == Float64
    @test size(grid) == (200,100,50)
    @test npoints(grid) == 200*100*50
    @test coordinates(grid, 1) == [0.,0.,0.]
    @test origin(grid) == (0.,0.,0.)
    @test spacing(grid) == (1.,1.,1.)

    grid = RegularGrid((-1.,-1.), (1.,1.), dims=(200,100))
    @test ndims(grid) == 2
    @test coordtype(grid) == Float64
    @test size(grid) == (200,100)
    @test npoints(grid) == 200*100
    @test coordinates(grid, 1) == [-1.,-1.]
    @test coordinates(grid, 200*100) == [1.,1.]
    @test origin(grid) == (-1.,-1.)

    grid = RegularGrid{Float64}(100,200)
    @test sprint(show, grid) == "100×200 RegularGrid{Float64,2}"
    @test sprint(show, MIME"text/plain"(), grid) == "RegularGrid{Float64,2}\n  dimensions: (100, 200)\n  origin:     (0.0, 0.0)\n  spacing:    (1.0, 1.0)"

    if visualtests
      gr(size=(800,800))
      @plottest plot(RegularGrid{Float64}(10)) joinpath(datadir,"grid1D.png") !istravis
      @plottest plot(RegularGrid{Float64}(10,20)) joinpath(datadir,"grid2D.png") !istravis
      @plottest plot(RegularGrid{Float64}(10,20,30)) joinpath(datadir,"grid3D.png") !istravis
      @plottest plot(RegularGrid{Float64}(10),[1,2,3,4,5,5,4,3,2,1]) joinpath(datadir,"grid1D-data.png") !istravis
      @plottest plot(RegularGrid{Float64}(10,10),1:100) joinpath(datadir,"grid2D-data.png") !istravis
      # @plottest plot(RegularGrid{Float64}(10,10,10),zeros(1000)) joinpath(datadir,"grid3D-data.png") !istravis
    end
  end

  @testset "StructuredGrid" begin
    nx, ny, nz = 20, 10, 10

    X = [1.,5.,6.,10.]
    g1 = StructuredGrid(X)
    @test npoints(g1) == 4
    @test coordinates(g1) == [1. 5. 6. 10.]

    @test sprint(show, g1) == "4 StructuredGrid{Float64,1}"
    @test sprint(show, MIME"text/plain"(), g1) == "4 StructuredGrid{Float64,1}\n 1.0  5.0  6.0  10.0"

    X  = [x for x in range(0,10,length=nx), j in 1:ny]
    Y  = sin.(X) .+ [0.5j for i in 1:nx, j in 1:ny]
    g2 = StructuredGrid(X, Y)
    @test npoints(g2) == nx*ny
    @test size(g2) == (nx, ny)
    @test size(coordinates(g2)) == (2, nx*ny)

    @test sprint(show, g2) == "20×10 StructuredGrid{Float64,2}"

    X = [x for x in range(0,10,length=nx), j in 1:ny, k in 1:nz]
    Y = sin.(X) .+ [0.5j for i in 1:nx, j in 1:ny, k in 1:nz]
    Z = [1.0(k-1) for i in 1:nx, j in 1:ny, k in 1:nz]
    g3 = StructuredGrid(X, Y, Z)
    @test npoints(g3) == nx*ny*nz
    @test size(g3) == (nx, ny, nz)
    @test size(coordinates(g3)) == (3, nx*ny*nz)

    @test sprint(show, g3) == "20×10×10 StructuredGrid{Float64,3}"

    if visualtests
      gr(size=(800,800))
      @plottest plot(g1) joinpath(datadir,"sgrid1D.png") !istravis
      @plottest plot(g2) joinpath(datadir,"sgrid2D.png") !istravis
      @plottest plot(g3,camera=(30,60)) joinpath(datadir,"sgrid3D.png") !istravis
      @plottest plot(g1,[1.,2.,2.,1]) joinpath(datadir,"sgrid1D-data.png") !istravis
      @plottest plot(g2,zeros(nx*ny)) joinpath(datadir,"sgrid2D-data.png") !istravis
      @plottest plot(g3,zeros(nx*ny*nz)) joinpath(datadir,"sgrid3D-data.png") !istravis
    end
  end
end


text:
function [fig] = PlotDG(c,CG,Trans,Topo,Param,fig,Slice)
OrdPoly=CG.OrdPoly;
figure(fig)
clf(fig)
hold on;
if strcmp(Param.Grid.Form,'Sphere') && Param.Flat
  dTol=2*pi/max(Param.nPanel-1,1);
end
switch Slice.Type
  case 'XY'
    level=Slice.iz;
    for iF=1:Param.Grid.NumFaces
      X=zeros(size(Param.Grid.Faces(iF).N,2),3);
      cLoc=reshape(c(:,:,iF),OrdPoly+1,OrdPoly+1);
      if strcmp(Param.Grid.Type,'Quad')
        dd=2/OrdPoly;
        ksi0=-1;
        for iRef=1:OrdPoly
          ksi1=ksi0+dd;
          eta0=-1;
          for jRef=1:OrdPoly
            eta1=eta0+dd;
            cc=0;
            for j=1:OrdPoly+1
              for i=1:OrdPoly+1
                cc=cc+Lagrange(0.5*(ksi0+ksi1),CG.xw,i)*Lagrange(0.5*(eta0+eta1),CG.xw,j)...
                  *cLoc(i,j);
              end
            end

            [X(1,:)]=Trans(ksi0,eta0,Param.Grid.z(level)...
              ,Param.Grid.Faces(iF),Topo,Param);
            [X(2,:)]=Trans(ksi1,eta0,Param.Grid.z(level)...
              ,Param.Grid.Faces(iF),Topo,Param);
            [X(3,:)]=Trans(ksi1,eta1,Param.Grid.z(level)...
              ,Param.Grid.Faces(iF),Topo,Param);
            [X(4,:)]=Trans(ksi0,eta1,Param.Grid.z(level)...
              ,Param.Grid.Faces(iF),Topo,Param);

            if strcmp(Param.Grid.Form,'Sphere') && Param.Flat
              [lam(1),theta(1)]=cart2sphere(X(1,1),X(1,2),X(1,3));
              [lam(2),theta(2)]=cart2sphere(X(2,1),X(2,2),X(2,3));
              [lam(3),theta(3)]=cart2sphere(X(3,1),X(3,2),X(3,3));
              [lam(4),theta(4)]=cart2sphere(X(4,1),X(4,2),X(4,3));
              lammin=min(lam);
              lammax=max(lam);
              if abs(lammin-lammax)>2*pi-dTol
                for i=1:4
                  if lam(i)<pi
                    lam(i)=lam(i)+2*pi;
                    if lam(i)>3*pi
                      lam(i)=lam(i)-2*pi;
                    end
                  end
                end
              end
              fill(lam,theta,real(cc));
            else
              fill3(X(:,1),X(:,2),X(:,3),real(cc));
            end

            eta0=eta1;
          end
          ksi0=ksi1;
        end
      else
      end
    end
  case 'XZ'
    y=Slice.y;
    for iF=1:Param.Grid.NumFaces
      for iz=1:Param.Grid.nz
        if y==Param.Grid.Faces(iF).Mid(2)
          cLoc=reshape(c(CG.Faces(iF).Glob,iz),OrdPoly+1,OrdPoly+1);
          dd=2/OrdPoly;
          eta0=-1;
          ksi0=-1;
          for iRef=1:OrdPoly
            ksi1=ksi0+dd;
            cc=0;
            for j=1:OrdPoly+1
              for i=1:OrdPoly+1
                cc=cc+Lagrange(0.5*(ksi0+ksi1),CG.xw,i)*Lagrange((eta0),CG.xw,j)...
                  *cLoc(i,j);
              end
            end
            [X(1,:)]=Trans(ksi0,eta0...
              ,Param.Grid.z(iz),Param.Grid.Faces(iF),Topo,Param);
            [X(2,:)]=Trans(ksi1,eta0...
              ,Param.Grid.z(iz),Param.Grid.Faces(iF),Topo,Param);
            [X(3,:)]=Trans(ksi1,eta0...
              ,Param.Grid.z(iz+1),Param.Grid.Faces(iF),Topo,Param);
            [X(4,:)]=Trans(ksi0,eta0...
              ,Param.Grid.z(iz+1),Param.Grid.Faces(iF),Topo,Param);
            fill3(X(:,1),X(:,2),X(:,3),cc);
            ksi0=ksi1;
          end
        end
      end
    end
  case 'YZ'
    x=Slice.x;
    for iF=1:Param.Grid.NumFaces
      for iz=1:Param.Grid.nz
        if x==Param.Grid.Faces(iF).Mid(1)
          cLoc=reshape(c(CG.Faces(iF).Glob,iz),OrdPoly+1,OrdPoly+1);
          dd=2/OrdPoly;
          eta0=-1;
          ksi0=-1;
          for jRef=1:OrdPoly
            eta1=eta0+dd;
            cc=0;
            for j=1:OrdPoly+1
              for i=1:OrdPoly+1
                cc=cc+Lagrange(ksi0,CG.xw,i)*Lagrange(0.5*(eta0+eta1),CG.xw,j)...
                  *cLoc(i,j);
              end
            end
            [X(1,:)]=Trans(ksi0,eta0...
              ,Param.Grid.z(iz),Param.Grid.Faces(iF),Topo,Param);
            [X(2,:)]=Trans(ksi0,eta1...
              ,Param.Grid.z(iz),Param.Grid.Faces(iF),Topo,Param);
            [X(3,:)]=Trans(ksi0,eta1...
              ,Param.Grid.z(iz+1),Param.Grid.Faces(iF),Topo,Param);
            [X(4,:)]=Trans(ksi0,eta0...
              ,Param.Grid.z(iz+1),Param.Grid.Faces(iF),Topo,Param);
            fill3(X(:,1),X(:,2),X(:,3),cc);
            eta0=eta1;
          end
        end
      end
    end
end
fig=fig+1;
hold off
end


text:
module DiffEqUncertainty

using DiffEqBase, Statistics, Distributions, Reexport
@reexport using Quadrature
include("probints.jl")
include("koopman.jl")

export ProbIntsUncertainty,AdaptiveProbIntsUncertainty
export expectation, centralmoment, Koopman, MonteCarlo

end


text:
# Autogenerated wrapper script for cmocka_jll for x86_64-apple-darwin
export libcmocka

JLLWrappers.@generate_wrapper_header("cmocka")
JLLWrappers.@declare_library_product(libcmocka, "@rpath/libcmocka.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcmocka,
        "lib/libcmocka.0.7.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#************************************************************************
# Global definitions
#************************************************************************

# Key and button actions
const RELEASE                = 0
const PRESS                  = 1
const REPEAT                 = 2

# Unknown key
const KEY_UNKNOWN            = -1

# Printable keys
const KEY_SPACE              = 32
const KEY_APOSTROPHE         = 39  # '
const KEY_COMMA              = 44  # ,
const KEY_MINUS              = 45  # -
const KEY_PERIOD             = 46  # .
const KEY_SLASH              = 47  # /
const KEY_0                  = 48
const KEY_1                  = 49
const KEY_2                  = 50
const KEY_3                  = 51
const KEY_4                  = 52
const KEY_5                  = 53
const KEY_6                  = 54
const KEY_7                  = 55
const KEY_8                  = 56
const KEY_9                  = 57
const KEY_SEMICOLON          = 59  # ;
const KEY_EQUAL              = 61  # =
const KEY_A                  = 65
const KEY_B                  = 66
const KEY_C                  = 67
const KEY_D                  = 68
const KEY_E                  = 69
const KEY_F                  = 70
const KEY_G                  = 71
const KEY_H                  = 72
const KEY_I                  = 73
const KEY_J                  = 74
const KEY_K                  = 75
const KEY_L                  = 76
const KEY_M                  = 77
const KEY_N                  = 78
const KEY_O                  = 79
const KEY_P                  = 80
const KEY_Q                  = 81
const KEY_R                  = 82
const KEY_S                  = 83
const KEY_T                  = 84
const KEY_U                  = 85
const KEY_V                  = 86
const KEY_W                  = 87
const KEY_X                  = 88
const KEY_Y                  = 89
const KEY_Z                  = 90
const KEY_LEFT_BRACKET       = 91  # [
const KEY_BACKSLASH          = 92  # \
const KEY_RIGHT_BRACKET      = 93  # ]
const KEY_GRAVE_ACCENT       = 96  # `
const KEY_WORLD_1            = 161 # non-US #1
const KEY_WORLD_2            = 162 # non-US #2

# Function keys
const KEY_ESCAPE             = 256
const KEY_ENTER              = 257
const KEY_TAB                = 258
const KEY_BACKSPACE          = 259
const KEY_INSERT             = 260
const KEY_DELETE             = 261
const KEY_RIGHT              = 262
const KEY_LEFT               = 263
const KEY_DOWN               = 264
const KEY_UP                 = 265
const KEY_PAGE_UP            = 266
const KEY_PAGE_DOWN          = 267
const KEY_HOME               = 268
const KEY_END                = 269
const KEY_CAPS_LOCK          = 280
const KEY_SCROLL_LOCK        = 281
const KEY_NUM_LOCK           = 282
const KEY_PRINT_SCREEN       = 283
const KEY_PAUSE              = 284
const KEY_F1                 = 290
const KEY_F2                 = 291
const KEY_F3                 = 292
const KEY_F4                 = 293
const KEY_F5                 = 294
const KEY_F6                 = 295
const KEY_F7                 = 296
const KEY_F8                 = 297
const KEY_F9                 = 298
const KEY_F10                = 299
const KEY_F11                = 300
const KEY_F12                = 301
const KEY_F13                = 302
const KEY_F14                = 303
const KEY_F15                = 304
const KEY_F16                = 305
const KEY_F17                = 306
const KEY_F18                = 307
const KEY_F19                = 308
const KEY_F20                = 309
const KEY_F21                = 310
const KEY_F22                = 311
const KEY_F23                = 312
const KEY_F24                = 313
const KEY_F25                = 314
const KEY_KP_0               = 320
const KEY_KP_1               = 321
const KEY_KP_2               = 322
const KEY_KP_3               = 323
const KEY_KP_4               = 324
const KEY_KP_5               = 325
const KEY_KP_6               = 326
const KEY_KP_7               = 327
const KEY_KP_8               = 328
const KEY_KP_9               = 329
const KEY_KP_DECIMAL         = 330
const KEY_KP_DIVIDE          = 331
const KEY_KP_MULTIPLY        = 332
const KEY_KP_SUBTRACT        = 333
const KEY_KP_ADD             = 334
const KEY_KP_ENTER           = 335
const KEY_KP_EQUAL           = 336
const KEY_LEFT_SHIFT         = 340
const KEY_LEFT_CONTROL       = 341
const KEY_LEFT_ALT           = 342
const KEY_LEFT_SUPER         = 343
const KEY_RIGHT_SHIFT        = 344
const KEY_RIGHT_CONTROL      = 345
const KEY_RIGHT_ALT          = 346
const KEY_RIGHT_SUPER        = 347
const KEY_MENU               = 348
const KEY_LAST               = KEY_MENU

# Modifier key bit flags
const MOD_SHIFT              = 0x0001
const MOD_CONTROL            = 0x0002
const MOD_ALT                = 0x0004
const MOD_SUPER              = 0x0008

# Mouse buttons
const MOUSE_BUTTON_1         = 0
const MOUSE_BUTTON_2         = 1
const MOUSE_BUTTON_3         = 2
const MOUSE_BUTTON_4         = 3
const MOUSE_BUTTON_5         = 4
const MOUSE_BUTTON_6         = 5
const MOUSE_BUTTON_7         = 6
const MOUSE_BUTTON_8         = 7
const MOUSE_BUTTON_LAST      = MOUSE_BUTTON_8
const MOUSE_BUTTON_LEFT      = MOUSE_BUTTON_1
const MOUSE_BUTTON_RIGHT     = MOUSE_BUTTON_2
const MOUSE_BUTTON_MIDDLE    = MOUSE_BUTTON_3

# Joysticks
const JOYSTICK_1             = 0
const JOYSTICK_2             = 1
const JOYSTICK_3             = 2
const JOYSTICK_4             = 3
const JOYSTICK_5             = 4
const JOYSTICK_6             = 5
const JOYSTICK_7             = 6
const JOYSTICK_8             = 7
const JOYSTICK_9             = 8
const JOYSTICK_10            = 9
const JOYSTICK_11            = 10
const JOYSTICK_12            = 11
const JOYSTICK_13            = 12
const JOYSTICK_14            = 13
const JOYSTICK_15            = 14
const JOYSTICK_16            = 15
const JOYSTICK_LAST          = JOYSTICK_16

# Error codes
const NOT_INITIALIZED        = 0x00010001  # GLFW has not been initialized.
const NO_CURRENT_CONTEXT     = 0x00010002  # No context is current for this thread.
const INVALID_ENUM           = 0x00010003  # One of the enum parameters for the function was given an invalid enum.
const INVALID_VALUE          = 0x00010004  # One of the parameters for the function was given an invalid value.
const OUT_OF_MEMORY          = 0x00010005  # A memory allocation failed.
const API_UNAVAILABLE        = 0x00010006  # GLFW could not find support for the requested client API on the system.
const VERSION_UNAVAILABLE    = 0x00010007  # The requested client API version is not available.
const PLATFORM_ERROR         = 0x00010008  # A platform-specific error occurred that does not match any of the more specific categories.
const FORMAT_UNAVAILABLE     = 0x00010009  # The clipboard did not contain data in the requested format.

const FOCUSED                = 0x00020001
const ICONIFIED              = 0x00020002
const RESIZABLE              = 0x00020003
const VISIBLE                = 0x00020004
const DECORATED              = 0x00020005
const AUTO_ICONIFY           = 0x00020006
const FLOATING               = 0x00020007

const RED_BITS               = 0x00021001
const GREEN_BITS             = 0x00021002
const BLUE_BITS              = 0x00021003
const ALPHA_BITS             = 0x00021004
const DEPTH_BITS             = 0x00021005
const STENCIL_BITS           = 0x00021006
const ACCUM_RED_BITS         = 0x00021007
const ACCUM_GREEN_BITS       = 0x00021008
const ACCUM_BLUE_BITS        = 0x00021009
const ACCUM_ALPHA_BITS       = 0x0002100A
const AUX_BUFFERS            = 0x0002100B
const STEREO                 = 0x0002100C
const SAMPLES                = 0x0002100D
const SRGB_CAPABLE           = 0x0002100E
const REFRESH_RATE           = 0x0002100F
const DOUBLEBUFFER           = 0x00021010

const CLIENT_API             = 0x00022001
const CONTEXT_VERSION_MAJOR  = 0x00022002
const CONTEXT_VERSION_MINOR  = 0x00022003
const CONTEXT_REVISION       = 0x00022004
const CONTEXT_ROBUSTNESS     = 0x00022005
const OPENGL_FORWARD_COMPAT  = 0x00022006
const OPENGL_DEBUG_CONTEXT   = 0x00022007
const OPENGL_PROFILE         = 0x00022008
const CONTEXT_RELEASE_BEHAVIOR = 0x00022009

const OPENGL_API             = 0x00030001
const OPENGL_ES_API          = 0x00030002

const NO_ROBUSTNESS          =          0
const NO_RESET_NOTIFICATION  = 0x00031001
const LOSE_CONTEXT_ON_RESET  = 0x00031002

const OPENGL_ANY_PROFILE     =          0
const OPENGL_CORE_PROFILE    = 0x00032001
const OPENGL_COMPAT_PROFILE  = 0x00032002

const CURSOR                 = 0x00033001
const STICKY_KEYS            = 0x00033002
const STICKY_MOUSE_BUTTONS   = 0x00033003

const CURSOR_NORMAL          = 0x00034001
const CURSOR_HIDDEN          = 0x00034002
const CURSOR_DISABLED        = 0x00034003

const ANY_RELEASE_BEHAVIOR   =          0
const RELEASE_BEHAVIOR_FLUSH = 0x00035001
const RELEASE_BEHAVIOR_NONE  = 0x00035002

# Standard cursor shapes
const ARROW_CURSOR           = 0x00036001
const IBEAM_CURSOR           = 0x00036002
const CROSSHAIR_CURSOR       = 0x00036003
const HAND_CURSOR            = 0x00036004
const HRESIZE_CURSOR         = 0x00036005
const VRESIZE_CURSOR         = 0x00036006

const CONNECTED              = 0x00040001
const DISCONNECTED           = 0x00040002

const DONT_CARE              = -1

#************************************************************************
# GLFW API types
#************************************************************************

immutable Monitor
	handle::Ptr{Void}
end
Base.show(io::IO, m::Monitor) = write(io, "Monitor($(m.handle == C_NULL ? m.handle : GetMonitorName(m)))")

typealias WindowHandle Ptr{Void}
type Window
	handle::WindowHandle
	callbacks::Vector{Function}
end
Window(handle::WindowHandle) = Window(handle, Vector{Function}(_window_callbacks_len))
import Base.==; ==(x::Window, y::Window) = (x.handle == y.handle)
Base.cconvert(::Type{WindowHandle}, window::Window) = window.handle
Base.cconvert(::Type{Window}, handle::WindowHandle) = ccall( (:glfwGetWindowUserPointer, lib), Ref{Window}, (WindowHandle,), handle)
Base.hash(window::Window, h::UInt64) = hash(window.handle, h)

immutable Cursor
	handle::Ptr{Void}
end

immutable VidMode
	width::Cint         # The width, in screen coordinates, of the video mode.
	height::Cint        # The height, in screen coordinates, of the video mode.
	redbits::Cint       # The bit depth of the red channel of the video mode.
	greenbits::Cint     # The bit depth of the green channel of the video mode.
	bluebits::Cint      # The bit depth of the blue channel of the video mode.
	refreshrate::Cint   # The refresh rate, in Hz, of the video mode.
end

#************************************************************************
# GLFW API functions
#************************************************************************

# Initialization and version information
Init() = Bool(ccall( (:glfwInit, lib), Cint, ())) || error("initialization failed")
Terminate() = ccall( (:glfwTerminate, lib), Void, ())
GetVersionString() = unsafe_string(ccall( (:glfwGetVersionString, lib), Cstring, ()))

# Error handling
@callback Error(code::Cint, description::Cstring) -> (code, unsafe_string(description))

# Monitor handling
function GetMonitors()
	count = Ref{Cint}()
	ptr = ccall( (:glfwGetMonitors, lib), Ptr{Monitor}, (Ref{Cint},), count)
	unsafe_wrap(Array, ptr, count[])
end

GetPrimaryMonitor() = ccall( (:glfwGetPrimaryMonitor, lib), Monitor, ())

function GetMonitorPos(monitor::Monitor)
	xpos, ypos = Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetMonitorPos, lib), Void, (Monitor, Ref{Cint}, Ref{Cint}), monitor, xpos, ypos)
	(xpos[], ypos[])
end

function GetMonitorPhysicalSize(monitor::Monitor)
	width, height = Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetMonitorPhysicalSize, lib), Void, (Monitor, Ref{Cint}, Ref{Cint}), monitor, width, height)
	(width[], height[])
end

GetMonitorName(monitor::Monitor) = unsafe_string(ccall( (:glfwGetMonitorName, lib), Cstring, (Monitor,), monitor))
@callback Monitor(monitor::Monitor, event::Cint)

function GetVideoModes(monitor::Monitor)
	count = Ref{Cint}()
	ptr = ccall( (:glfwGetVideoModes, lib), Ptr{VidMode}, (Monitor, Ref{Cint}), monitor, count)
	unsafe_wrap(Array, ptr, count[])
end

GetVideoMode(monitor::Monitor) = unsafe_load(ccall( (:glfwGetVideoMode, lib), Ptr{VidMode}, (Monitor,), monitor))
SetGamma(monitor::Monitor, gamma::Real) = ccall( (:glfwSetGamma, lib), Void, (Monitor, Cfloat), monitor, gamma)

# Window handling
DefaultWindowHints() = ccall( (:glfwDefaultWindowHints, lib), Void, ())
WindowHint(target::Integer, hint::Integer) = ccall( (:glfwWindowHint, lib), Void, (Cint, Cint), target, hint)

function CreateWindow(width::Integer, height::Integer, title::AbstractString, monitor::Monitor=Monitor(C_NULL), share::Window=Window(C_NULL))
	handle = ccall( (:glfwCreateWindow, lib), WindowHandle, (Cint, Cint, Cstring, Monitor, WindowHandle), width, height, title, monitor, share)
	window = Window(handle)
	ccall( (:glfwSetWindowUserPointer, lib), Void, (WindowHandle, Ref{Window}), window, window)
	window
end

DestroyWindow(window::Window) = (yield(); ccall( (:glfwDestroyWindow, lib), Void, (WindowHandle,), window))
WindowShouldClose(window::Window) = ccall( (:glfwWindowShouldClose, lib), Cint, (WindowHandle,), window) != 0
SetWindowShouldClose(window::Window, value::Bool) = ccall( (:glfwSetWindowShouldClose, lib), Void, (WindowHandle, Cint), window, value)
SetWindowTitle(window::Window, title::AbstractString) = ccall( (:glfwSetWindowTitle, lib), Void, (WindowHandle, Cstring), window, title)

function GetWindowPos(window::Window)
	xpos, ypos = Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetWindowPos, lib), Void, (WindowHandle, Ref{Cint}, Ref{Cint}), window, xpos, ypos)
	(xpos[], ypos[])
end

SetWindowPos(window::Window, xpos::Integer, ypos::Integer) = ccall( (:glfwSetWindowPos, lib), Void, (WindowHandle, Cint, Cint), window, xpos, ypos)

function GetWindowSize(window::Window)
	width, height = Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetWindowSize, lib), Void, (WindowHandle, Ref{Cint}, Ref{Cint}), window, width, height)
	(width[], height[])
end

SetWindowSize(window::Window, width::Integer, height::Integer) = ccall( (:glfwSetWindowSize, lib), Void, (WindowHandle, Cint, Cint), window, width, height)

function GetFramebufferSize(window::Window)
	width, height = Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetFramebufferSize, lib), Void, (WindowHandle, Ref{Cint}, Ref{Cint}), window, width, height)
	(width[], height[])
end

function GetWindowFrameSize(window::Window)
	left, top, right, bottom = Ref{Cint}(), Ref{Cint}(), Ref{Cint}(), Ref{Cint}()
	ccall( (:glfwGetWindowFrameSize, lib), Void, (WindowHandle, Ref{Cint}, Ref{Cint}, Ref{Cint}, Ref{Cint}), window, left, top, right, bottom)
	(left[], top[], right[], bottom[])
end

IconifyWindow(window::Window) = ccall( (:glfwIconifyWindow, lib), Void, (WindowHandle,), window)
RestoreWindow(window::Window) = ccall( (:glfwRestoreWindow, lib), Void, (WindowHandle,), window)
ShowWindow(window::Window) = ccall( (:glfwShowWindow, lib), Void, (WindowHandle,), window)
HideWindow(window::Window) = ccall( (:glfwHideWindow, lib), Void, (WindowHandle,), window)
GetWindowMonitor(window::Window) = ccall( (:glfwGetWindowMonitor, lib), Monitor, (WindowHandle,), window)
GetWindowAttrib(window::Window, attrib::Integer) = ccall( (:glfwGetWindowAttrib, lib), Cint, (WindowHandle, Cint), window, attrib)
@callback WindowPos(window::Window, xpos::Cint, ypos::Cint)
@callback WindowSize(window::Window, width::Cint, height::Cint)
@callback WindowClose(window::Window)
@callback WindowRefresh(window::Window)
@callback WindowFocus(window::Window, focused::Cint) -> (window, Bool(focused))
@callback WindowIconify(window::Window, iconified::Cint) -> (window, Bool(iconified))
@callback FramebufferSize(window::Window, width::Cint, height::Cint)
PollEvents() = ccall( (:glfwPollEvents, lib), Void, ())
WaitEvents() = ccall( (:glfwWaitEvents, lib), Void, ())
PostEmptyEvent() = ccall( (:glfwPostEmptyEvent, lib), Void, ())

# Input handling
function GetInputMode(window::Window, mode::Integer)
	value = ccall( (:glfwGetInputMode, lib), Cint, (WindowHandle, Cint), window, mode)
	if mode in (STICKY_KEYS, STICKY_MOUSE_BUTTONS)
		value = Bool(value)
	end
	value
end

SetInputMode(window::Window, mode::Integer, value::Integer) = ccall( (:glfwSetInputMode, lib), Void, (WindowHandle, Cint, Cint), window, mode, value)
GetKey(window::Window, key::Integer) = Bool(ccall( (:glfwGetKey, lib), Cint, (WindowHandle, Cint), window, key))
GetMouseButton(window::Window, button::Integer) = Bool(ccall( (:glfwGetMouseButton, lib), Cint, (WindowHandle, Cint), window, button))

function GetCursorPos(window::Window)
	xpos, ypos = Ref{Cdouble}(), Ref{Cdouble}()
	ccall( (:glfwGetCursorPos, lib), Void, (WindowHandle, Ref{Cdouble}, Ref{Cdouble}), window, xpos, ypos)
	(xpos[], ypos[])
end

SetCursorPos(window::Window, xpos::Real, ypos::Real) = ccall( (:glfwSetCursorPos, lib), Void, (WindowHandle, Cdouble, Cdouble), window, xpos, ypos)
CreateStandardCursor(shape::Integer) = ccall( (:glfwCreateStandardCursor, lib), Cursor, (Cint,), shape)
DestroyCursor(cursor::Cursor) = ccall( (:glfwDestroyCursor, lib), Void, (Cursor,), cursor)
SetCursor(window::Window, cursor::Cursor) = ccall( (:glfwSetCursor, lib), Void, (WindowHandle, Cursor), window, cursor)
SetCursor(window::Window, ::Void) = SetCursor(window, Cursor(C_NULL))
@callback Key(window::Window, key::Cint, scancode::Cint, action::Cint, mods::Cint)
@callback Char(window::Window, codepoint::Cuint) -> (window, convert(Char, codepoint))
@callback CharMods(window::Window, codepoint::Cuint, mods::Cint) -> (window, convert(Char, codepoint), mods)
@callback MouseButton(window::Window, button::Cint, action::Cint, mods::Cint)
@callback CursorPos(window::Window, xpos::Cdouble, ypos::Cdouble)
@callback CursorEnter(window::Window, entered::Cint) -> (window, Bool(entered))
@callback Scroll(window::Window, xoffset::Cdouble, yoffset::Cdouble)
@callback Drop(window::Window, count::Cint, paths::Ptr{Cstring}) -> (window, map(unsafe_string, unsafe_wrap(Array, paths, count)))
JoystickPresent(joy::Integer) = Bool(ccall( (:glfwJoystickPresent, lib), Cint, (Cint,), joy))

function GetJoystickAxes(joy::Integer)
	count = Ref{Cint}()
	ptr = ccall( (:glfwGetJoystickAxes, lib), Ptr{Cfloat}, (Cint, Ref{Cint}), joy, count)
	unsafe_wrap(Array, ptr, count[])
end

function GetJoystickButtons(joy::Integer)
	count = Ref{Cint}()
	ptr = ccall( (:glfwGetJoystickButtons, lib), Ptr{Cuchar}, (Cint, Ref{Cint}), joy, count)
	unsafe_wrap(Array, ptr, count[])
end

GetJoystickName(joy::Integer) = unsafe_string(ccall( (:glfwGetJoystickName, lib), Cstring, (Cint,), joy))

# Context handling
MakeContextCurrent(window::Window) = ccall( (:glfwMakeContextCurrent, lib), Void, (WindowHandle,), window)
GetCurrentContext() = Base.cconvert(Window, ccall( (:glfwGetCurrentContext, lib), WindowHandle, ()))
SwapBuffers(window::Window) = ccall( (:glfwSwapBuffers, lib), Void, (WindowHandle,), window)
SwapInterval(interval::Integer) = ccall( (:glfwSwapInterval, lib), Void, (Cint,), interval)
ExtensionSupported(extension::AbstractString) = Bool(ccall( (:glfwExtensionSupported, lib), Cint, (Cstring,), extension))
GetProcAddress(procname::AbstractString) = ccall((:glfwGetProcAddress, lib), Ptr{Void}, (Cstring,), procname)


text:
using CUDex

x = DexArray(rand(6,5,4,3))
w = DexArray(rand(3,3,4,2))

t1 = conv_f(x, w)
t2 = relu(t1)
t3 = pool_f(t2)


text:
init_wrapper_func(spec::SpecFunc) =
    Dict(:category => :function, :name => nc_convert(SnakeCaseLower, remove_vk_prefix(spec.name)), :short => false)
init_wrapper_func(spec::Spec) = Dict(:category => :function, :name => remove_vk_prefix(spec.name), :short => false)

arg_decl(x::Spec) = :($(wrap_identifier(x))::$(signature_type(nice_julian_type(x))))

function kwarg_decl(x::Spec)
    val = @match x begin
        GuardBy(is_default_count) => default_count(x)
        _ => default(x)
    end
    Expr(:kw, wrap_identifier(x), val)
end

drop_arg(x::Spec) = is_length(x) && !is_specific_count(x) && !is_default_count(x) || is_pointer_start(x) || x.type == :(Ptr{Ptr{Cvoid}})

"""
Function pointer arguments for a handle.
Includes one `fun_ptr_create` for the constructor (if applicable),
and one `fun_ptr_destroy` for the destructor (if applicable).
"""
function func_ptr_args(spec::SpecHandle)
    args = Expr[]
    spec ∈ spec_create_funcs.handle && push!(args, :(fun_ptr_create::FunctionPtr))
    destructor(spec) ≠ :identity && push!(args, :(fun_ptr_destroy::FunctionPtr))
    args
end

"""
Function pointer arguments for a function.
Takes the function pointers arguments of the underlying handle if it is a Vulkan constructor,
or a unique `fun_ptr` if that's just a normal Vulkan function.
"""
function func_ptr_args(spec::SpecFunc)
    if spec.type ∈ [CREATE, ALLOCATE]
        func_ptr_args(create_func(spec).handle)
    else
        [:(fun_ptr::FunctionPtr)]
    end
end

"""
Corresponding pointer argument for a Vulkan function.
"""
func_ptrs(spec::Spec) = name.(func_ptr_args(spec))

function add_func_args!(p::Dict, spec, params; with_func_ptr = false)
    params = filter(!drop_arg, params)
    arg_filter = if spec.type ∈ [DESTROY, FREE]
        destroyed_type = destroy_func(spec).handle.name
        x -> !is_optional(x) || x.type == destroyed_type
    else
        !is_optional
    end

    p[:args] = map(arg_decl, filter(arg_filter, params))
    p[:kwargs] = map(kwarg_decl, filter(!arg_filter, params))

    with_func_ptr && append!(p[:args], func_ptr_args(spec))
end


text:
##############

Base.@kwdef mutable struct Markers
    """
        Data structure for marker fields
    """
    x         ::  Union{Array{Float64,1}, Missing} = missing
    y         ::  Union{Array{Float64,1}, Missing} = missing
    T         ::  Union{Array{Float64,1}, Missing} = missing
    phase     ::  Union{Array{Float64,1}, Missing} = missing
    cellx     ::  Union{Array{Int64,1}, Missing} = missing
    celly     ::  Union{Array{Int64,1}, Missing} = missing
    nmark     ::  Union{Int64, Missing} = missing
    nmark_ini ::  Union{Int64, Missing} = missing
    nmark_max ::  Union{Int64, Missing} = missing
    nmx       ::  Union{Int64, Missing} = missing
    nmy       ::  Union{Int64, Missing} = missing
end
export Markers

##############

Base.@kwdef mutable struct Thermal_BC
    type_W  ::Union{Int64, Missing} = 1 # Default Dirichlet
    type_E  ::Union{Int64, Missing} = 1 # Default Dirichlet
    type_S  ::Union{Int64, Missing} = 1 # Default Dirichlet
    type_N  ::Union{Int64, Missing} = 1 # Default Dirichlet
    Dir_W   ::Union{Vector{Float64}, Missing} = missing
    Dir_E   ::Union{Vector{Float64}, Missing} = missing
    Dir_S   ::Union{Vector{Float64}, Missing} = missing
    Dir_N   ::Union{Vector{Float64}, Missing} = missing
    Pinned  ::Union{Matrix{Int64}, Missing}   = missing
end
export Thermal_BC

##############

Base.@kwdef mutable struct BoundaryConditions
    T   :: Union{Thermal_BC, Missing} = Thermal_BC()
    Vx  :: Union{Thermal_BC, Missing} = Thermal_BC()
    Vy  :: Union{Thermal_BC, Missing} = Thermal_BC()
    periodix  ::Union{Int64, Missing} = 0
end
export BoundaryConditions

##############

Base.@kwdef mutable struct Tensor2D
    xx   :: Union{Matrix{Float64}, Missing} = missing
    yy   :: Union{Matrix{Float64}, Missing} = missing
    zz   :: Union{Matrix{Float64}, Missing} = missing
    xy   :: Union{Matrix{Float64}, Missing} = missing
    xy_c :: Union{Matrix{Float64}, Missing} = missing
    II   :: Union{Matrix{Float64}, Missing} = missing
    div  :: Union{Matrix{Float64}, Missing} = missing
end
export Tensor2D

##############

Base.@kwdef mutable struct ModelDomain
    xc   :: Union{LinRange{Float64}, Missing} = missing
    yc   :: Union{LinRange{Float64}, Missing} = missing
    xv   :: Union{LinRange{Float64}, Missing} = missing
    yv   :: Union{LinRange{Float64}, Missing} = missing
    xce  :: Union{LinRange{Float64}, Missing} = missing
    yce  :: Union{LinRange{Float64}, Missing} = missing
    xmin :: Union{Float64, Missing} = missing
    xmax :: Union{Float64, Missing} = missing
    ymin :: Union{Float64, Missing} = missing
    ymax :: Union{Float64, Missing} = missing
    dx   :: Union{Float64, Missing} = missing
    dy   :: Union{Float64, Missing} = missing
    ncx  :: Union{Int64, Missing}   = missing
    ncy  :: Union{Int64, Missing}   = missing
    L    :: Union{Float64, Missing} = missing
    L0   :: Union{Float64, Missing} = missing
end
export ModelDomain

##############

Base.@kwdef mutable struct MaterialParameters
    nphase       :: Union{Int64, Missing}   = missing
    # Viscosity
    n    :: Union{Vector{Float64}, Missing} = missing
    eta0 :: Union{Vector{Float64}, Missing} = missing
    Tref :: Union{Vector{Float64}, Missing} = missing
    # Elasticity
    K    :: Union{Vector{Float64}, Missing} = missing
    G    :: Union{Vector{Float64}, Missing} = missing
    # Phase field damage
    # lDam :: Union{Vector{Float64}, Missing} = missing
    # gDam :: Union{Vector{Float64}, Missing} = missing
    # eDam :: Union{Vector{Float64}, Missing} = missing
    lDam :: Union{Float64, Missing} = missing
    gDam :: Union{Float64, Missing} = missing
    eDam :: Union{Float64, Missing} = missing
end
export MaterialParameters

##############

Base.@kwdef mutable struct ModelParameters
    # User interface
    user          = zeros(Float64, 10)
    # Switches
    BC_type       :: Union{Int64, Missing}     = missing    # 1: Pure shear / 2: Simple shear / 3: Simple shear periodic
    comp          :: Union{Int64, Missing}     = missing    # Compressible
    # Time integration
    advection     :: Union{Int64, Missing}     = missing
    ALE           :: Union{Int64, Missing}     = missing    # Deform box
    Ebg           :: Union{Float64, Missing}   = missing
    nt            :: Union{Int64, Missing}     = missing
    dt            :: Union{Float64, Missing}   = missing
    Courant       :: Union{Float64, Missing}   = missing
    dt_var        :: Union{Int64, Missing}     = missing
    t             :: Union{Float64, Missing}   = missing
    rkw           = 1.0/6.0*[1.0 2.0 2.0 1.0] # RK4 weights for averaging
    rkv           = 1.0/2.0*[1.0 1.0 2.0 2.0] # RK4 weights for time stepping
    # Linear solver
    solver        :: Union{Int64, Missing}     = missing    # 0: coupled --- 1: decoupled
    gamma         :: Union{Float64, Missing}   = missing    # Penalty factor
    Dir_scale     :: Union{Float64, Missing}   = missing    # Dirichlet scaling factor
    # Non-linear solver
    niter_nl      :: Union{Int64, Missing}     = missing    # max. number of non-linear iterations
    tol_nl        :: Union{Float64, Missing}   = missing    # non-linear tolerance
    JFNK          :: Union{Int64, Missing}     = missing    # Jacobian-Free Newton_Krylov
    # Visualisation
    show_figs     :: Union{Int64, Missing}     = missing    # activates visualisation...
    nout          :: Union{Int64, Missing}     = missing    # ... every nout
end
export ModelParameters

##############

Base.@kwdef mutable struct Fields2D
    etac      :: Union{Matrix{Float64}, Missing} = missing
    etav      :: Union{Matrix{Float64}, Missing} = missing
    Kc        :: Union{Matrix{Float64}, Missing} = missing
    Gc        :: Union{Matrix{Float64}, Missing} = missing
    Gv        :: Union{Matrix{Float64}, Missing} = missing
    We        :: Union{Matrix{Float64}, Missing} = missing
    We0       :: Union{Matrix{Float64}, Missing} = missing
    rhoc      :: Union{Matrix{Float64}, Missing} = missing
    Pc0       :: Union{Matrix{Float64}, Missing} = missing
    Pc        :: Union{Matrix{Float64}, Missing} = missing
    Vx        :: Union{Matrix{Float64}, Missing} = missing
    Vy        :: Union{Matrix{Float64}, Missing} = missing
    Vxc       :: Union{Matrix{Float64}, Missing} = missing
    Vyc       :: Union{Matrix{Float64}, Missing} = missing
    div       :: Union{Matrix{Float64}, Missing} = missing
    Fx        :: Union{Matrix{Float64}, Missing} = missing
    Fy        :: Union{Matrix{Float64}, Missing} = missing
    Fp        :: Union{Matrix{Float64}, Missing} = missing
    FDam      :: Union{Matrix{Float64}, Missing} = missing
    Damc      :: Union{Matrix{Float64}, Missing} = missing
    Damv      :: Union{Matrix{Float64}, Missing} = missing
    phiDam    :: Union{Matrix{Float64}, Missing} = missing
    phiDam0   :: Union{Matrix{Float64}, Missing} = missing
    GDam      :: Union{Matrix{Float64}, Missing} = missing
    lDam      :: Union{Matrix{Float64}, Missing} = missing
    #####
    NumVx     :: Union{Matrix{Int64}, Missing} = missing
    NumVy     :: Union{Matrix{Int64}, Missing} = missing
    NumP      :: Union{Matrix{Int64}, Missing} = missing
    #####
    mpc           :: Union{Matrix{Float64}, Missing} = missing
    mpc_th        :: Union{Array{Float64, 3}, Missing} = missing
    phase_perc    :: Union{Array{Float64, 3}, Missing} = missing
    phase_perc_th :: Union{Array{Float64, 4}, Missing} = missing
end
export Fields2D

##############

Base.@kwdef mutable struct CharScales
    L      :: Union{Float64, Missing}   = missing
    t      :: Union{Float64, Missing}   = missing
    S      :: Union{Float64, Missing}   = missing
    T      :: Union{Float64, Missing}   = missing
end
export CharScales


text:
export Forest

mutable struct Forest{T<:AbstractFloat} <: Environment
  worldAABBmin::Vector{T}
  worldAABBmax::Vector{T}
  keepin_zones::Vector
  keepout_zones::Vector
  obstacle_set::Vector
  obstacles_set_render::Vector
end

function Forest{T}() where T
  vars = matread(joinpath(Pkg.dir("GuSTO"), "src", "environment","forest.mat"))

  keepout_zones = Vector{HyperRectangle}(0)
  push!(keepout_zones, HyperRectangle(Vec3f0(-20.,-20.,28.), Vec3f0(160.,160.,20)))
  push!(keepout_zones, HyperRectangle(Vec3f0(-20.,-20.,-20.), Vec3f0(160.,160.,20)))
  push!(keepout_zones, HyperRectangle(Vec3f0(-20.,-20.,-20.), Vec3f0(20.,160.,68)))
  push!(keepout_zones, HyperRectangle(Vec3f0(120.,-20.,-20.), Vec3f0(20.,160.,68)))
  push!(keepout_zones, HyperRectangle(Vec3f0(-20.,-20.,-20.), Vec3f0(160.,20.,68)))
  push!(keepout_zones, HyperRectangle(Vec3f0(-20.,120.,-20.), Vec3f0(160.,20.,68)))

  obstacle_set = Vector{HyperRectangle}(0)
  for idx in 1:2:size(vars["obstacles_infl"],1)
    c1,c2 = vars["obstacles_infl"][idx,:], vars["obstacles_infl"][idx+1,:]
    push!(obstacle_set,
      HyperRectangle(Vec3f0(c1),Vec3f0(c2-c1)))
  end

  obstacles_set_render = Vector{HyperRectangle}(0)
  keys = ["tree_obstacles", "tower_obstacles"]
  for key in keys
    for idx in 1:2:size(vars[key],1)
      c1,c2 = vars[key][idx,:], vars[key][idx+1,:]
      push!(obstacles_set_render,
        HyperRectangle(Vec3f0(c1),Vec3f0(c2-c1)))
    end
  end

  worldAABBmin = Inf*ones(T,3)
  worldAABBmax = -Inf*ones(T,3)
  for zone in keepout_zones
    corner1,corner2 = zone.origin, zone.origin+zone.widths
    zone_min = [min(corner1[i], corner2[i]) for i in 1:3]
    worldAABBmin[worldAABBmin .> zone_min] = zone_min[worldAABBmin .> zone_min]
    zone_max = [max(corner1[i], corner2[i]) for i in 1:3]
    worldAABBmax[worldAABBmax .< zone_max] = zone_max[worldAABBmax .< zone_max]
  end

  keepin_zones = Vector{HyperRectangle}(0)
  widths = worldAABBmax-worldAABBmin
  margin = 0.1*widths
  push!(keepin_zones, HyperRectangle(Vec3f0(worldAABBmin-margin),Vec3f0(widths+margin)))

  return Forest{T}(worldAABBmin, worldAABBmax, keepin_zones,
    keepout_zones, obstacle_set, obstacles_set_render)
end
Forest(::Type{T} = Float64; kwargs...) where {T} = Forest{T}(; kwargs...)

# function plot_forest(f::Forest)
#   function plot_shape(rect::HyperRectangle,plot_color::String="blue")
#     lenX,lenY,lenZ = rect.widths
#
#     (x1,y1,z1) = rect.origin
#     (x2,y2,z2) = rect.origin+rect.widths
#
#     xs = [x1; x2; x2; x1; x1]
#     ys = [y1; y1; y2; y2; y1]
#     zs = [z1; z1; z1; z1; z1]
#
#     zs2 = [z2; z2; z2; z2; z2]
#
#     xs3 = [x1; x1; x2; x2; x2; x2; x1; x1]
#     ys3 = [y2; y2; y2; y2; y1; y1; y1; y1]
#     zs3 = [z1; z2; z2; z1; z1; z2; z2; z1]
#
#     PyPlot.mesh(xs, ys, zs, color=plot_color)
#     PyPlot.mesh(xs, ys, zs2, color=plot_color)
#     PyPlot.mesh(xs3, ys3, zs3, color=plot_color)
#   end
#
#   PyPlot.figure()
#   plot_shape(f.keepin_zones[1],"green")
#
#   for z in f.keepout_zones
#     plot_shape(z, "red")
#   end
#   for z in f.obstacles_set_render
#     plot_shape(z, "black")
#   end
# end


text:
import ...BalanceLaws:
    tendency_from_slow_to_fast!,
    cummulate_fast_solution!,
    reconcile_from_fast_to_slow!

#using Printf

@inline function initialize_fast_state!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
    slow_dt,
    fast_time_rec,
    fast_steps,
)

    #- inverse ratio of additional fast time steps (for weighted average)
    #  --> do 1/add more time-steps and average from: 1 - 1/add up to: 1 + 1/add
    add = slow.add_fast_substeps

    #- set time-step and number of sub-steps we need
    #  will time-average fast over: fast_steps[1] , fast_steps[3]
    #  centered on fast_steps[2] which corresponds to advance in time of slow
    fast_dt = fast_time_rec[1]
    if add == 0
        steps = fast_dt > 0 ? ceil(Int, slow_dt / fast_dt) : 1
        fast_steps[1:3] = [1 1 1] * steps
    else
        steps = fast_dt > 0 ? ceil(Int, slow_dt / fast_dt / add) : 1
        fast_steps[2] = add * steps
        fast_steps[1] = (add - 1) * steps
        fast_steps[3] = (add + 1) * steps
    end
    fast_time_rec[1] = slow_dt / fast_steps[2]
    fast_time_rec[2] = 0.0
    # @printf("Update: frac_dt = %.1f , dt_fast = %.1f , nsubsteps= %i\n",
    #          slow_dt,fast_time_rec[1],fast_steps[3])
    # println(" fast_steps = ",fast_steps)

    dgFast.state_auxiliary.η_c .= -0
    dgFast.state_auxiliary.U_c .= (@SVector [-0, -0])'

    # preliminary test: no average
    Qfast.η .= dgFast.state_auxiliary.η_s
    Qfast.U .= dgFast.state_auxiliary.U_s

    return nothing
end

@inline function initialize_adjustment!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
)
    ## reset tendency adjustment before calling Baroclinic Model
    dgSlow.state_auxiliary.ΔGu .= 0

    return nothing
end

@inline function tendency_from_slow_to_fast!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
    dQslow2fast,
)
    FT = eltype(Qslow)

    # integrate the tendency
    tendency_dg = dgSlow.modeldata.tendency_dg
    tend = tendency_dg.balance_law
    grid = dgSlow.grid
    elems = grid.topology.elems
    update_auxiliary_state!(tendency_dg, tend, dQslow2fast, 0, elems)

    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(dgSlow)

    ## get top value (=integral over full depth) of ∫du
    nb_aux_tnd = number_states(tend, Auxiliary())
    data_tnd = reshape(
        tendency_dg.state_auxiliary.data,
        Nq^2,
        Nqk,
        nb_aux_tnd,
        nelemv,
        nelemh,
    )
    index_∫du = varsindex(vars_state(tend, Auxiliary(), FT), :∫du)
    flat_∫du = @view data_tnd[:, end, index_∫du, end, 1:nrealelemh]

    ## copy into Gᵁ of dgFast
    nb_aux_fst = number_states(fast, Auxiliary())
    data_fst = reshape(dgFast.state_auxiliary.data, Nq^2, nb_aux_fst, nelemh)
    index_Gᵁ = varsindex(vars_state(fast, Auxiliary(), FT), :Gᵁ)
    boxy_Gᵁ = @view data_fst[:, index_Gᵁ, 1:nrealelemh]
    boxy_Gᵁ .= flat_∫du

    ## scale by -1/H and copy back to ΔGu
    # note: since tendency_dg.state_auxiliary.∫du is not used after this, could be
    #   re-used to store a 3-D copy of "-Gu"
    nb_aux_slw = number_states(slow, Auxiliary())
    data_slw = reshape(
        dgSlow.state_auxiliary.data,
        Nq^2,
        Nqk,
        nb_aux_slw,
        nelemv,
        nelemh,
    )
    index_ΔGu = varsindex(vars_state(slow, Auxiliary(), FT), :ΔGu)
    boxy_ΔGu = @view data_slw[:, :, index_ΔGu, :, 1:nrealelemh]
    boxy_∫du = @view data_tnd[:, end:end, index_∫du, end:end, 1:nrealelemh]
    boxy_ΔGu .= -boxy_∫du / slow.problem.H

    return nothing
end

@inline function cummulate_fast_solution!(
    fast::BarotropicModel,
    dgFast,
    Qfast,
    fast_time,
    fast_dt,
    substep,
    fast_steps,
    fast_time_rec,
)
    #- might want to use some of the weighting factors: weights_η & weights_U
    #- should account for case where fast_dt < fast.param.dt

    # cumulate Fast solution:
    if substep >= fast_steps[1]
        dgFast.state_auxiliary.U_c .+= Qfast.U
        dgFast.state_auxiliary.η_c .+= Qfast.η
        fast_time_rec[2] += 1.0
    end

    # save mid-point solution to start from the next time-step
    if substep == fast_steps[2]
        dgFast.state_auxiliary.U_s .= Qfast.U
        dgFast.state_auxiliary.η_s .= Qfast.η
    end

    return nothing
end

@inline function reconcile_from_fast_to_slow!(
    slow::OceanModel,
    fast::BarotropicModel,
    dgSlow,
    dgFast,
    Qslow,
    Qfast,
    fast_time_rec,
)
    FT = eltype(Qslow)
    Nq, Nqk, _, _, nelemv, nelemh, nrealelemh, _ = basic_grid_info(dgSlow)
    grid = dgSlow.grid
    elems = grid.topology.elems

    # need to calculate int_u using integral kernels
    # u_slow := u_slow + (1/H) * (u_fast - \int_{-H}^{0} u_slow)

    ## get time weighted averaged out of cumulative arrays
    dgFast.state_auxiliary.U_c .*= 1 / fast_time_rec[2]
    dgFast.state_auxiliary.η_c .*= 1 / fast_time_rec[2]

    ## Compute: \int_{-H}^{0} u_slow)

    # integrate vertically horizontal velocity
    flowintegral_dg = dgSlow.modeldata.flowintegral_dg
    flowint = flowintegral_dg.balance_law
    update_auxiliary_state!(flowintegral_dg, flowint, Qslow, 0, elems)

    # get top value (=integral over full depth)
    nb_aux_flw = number_states(flowint, Auxiliary())
    data_flw = reshape(
        flowintegral_dg.state_auxiliary.data,
        Nq^2,
        Nqk,
        nb_aux_flw,
        nelemv,
        nelemh,
    )
    index_∫u = varsindex(vars_state(flowint, Auxiliary(), FT), :∫u)
    flat_∫u = @view data_flw[:, end, index_∫u, end, 1:nrealelemh]

    ## substract ∫u from U and divide by H

    # Δu is a place holder for 1/H * (Ū - ∫u)
    Δu = dgFast.state_auxiliary.Δu
    Δu .= dgFast.state_auxiliary.U_c

    nb_aux_fst = number_states(fast, Auxiliary())
    data_fst = reshape(dgFast.state_auxiliary.data, Nq^2, nb_aux_fst, nelemh)
    index_Δu = varsindex(vars_state(fast, Auxiliary(), FT), :Δu)
    boxy_Δu = @view data_fst[:, index_Δu, 1:nrealelemh]
    boxy_Δu .-= flat_∫u
    boxy_Δu ./= slow.problem.H

    ## apply the 2D correction to the 3D solution
    nb_cons_slw = number_states(slow, Prognostic())
    data_slw = reshape(Qslow.data, Nq^2, Nqk, nb_cons_slw, nelemv, nelemh)
    index_u = varsindex(vars_state(slow, Prognostic(), FT), :u)
    boxy_u = @view data_slw[:, :, index_u, :, 1:nrealelemh]
    boxy_u .+= reshape(boxy_Δu, Nq^2, 1, 2, 1, nrealelemh)

    ## save Eta from 3D model into η_diag (aux var of 2D model)
    ## and store difference between η from Barotropic Model and η_diag
    index_η = varsindex(vars_state(slow, Prognostic(), FT), :η)
    boxy_η_3D = @view data_slw[:, :, index_η, :, 1:nrealelemh]
    flat_η = @view data_slw[:, end, index_η, end, 1:nrealelemh]
    index_η_diag = varsindex(vars_state(fast, Auxiliary(), FT), :η_diag)
    boxy_η_diag = @view data_fst[:, index_η_diag, 1:nrealelemh]
    boxy_η_diag .= flat_η
    dgFast.state_auxiliary.Δη .=
        dgFast.state_auxiliary.η_c - dgFast.state_auxiliary.η_diag

    ## copy 2D model Eta over to 3D model
    index_η_c = varsindex(vars_state(fast, Auxiliary(), FT), :η_c)
    boxy_η_2D = @view data_fst[:, index_η_c, 1:nrealelemh]
    boxy_η_3D .= reshape(boxy_η_2D, Nq^2, 1, 1, 1, nrealelemh)

    return nothing
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA