text:
using Base
using Random
using Distributions

abstract type AbstractSyntheticDataSet end

# MAML Sinusoid
struct Sinusoid <: AbstractSyntheticDataSet
    p_A::UnivariateDistribution
    p_ϕ::UnivariateDistribution
    p_x::UnivariateDistribution
end
Sinusoid() = Sinusoid(Uniform(0.1, 5.0), Uniform(0.0, π), Uniform(-5.0, 5.0))

function Base.rand(p::Sinusoid, k::Int)
    A = rand(p.p_A)
    ϕ = rand(p.p_ϕ)
    X = rand(p.p_x, k)
    Y = @. A * sin(X + ϕ)
    return X, Y
end
function Base.rand(p::Sinusoid)
    A= rand(p.p_A)
    ϕ = rand(p.p_ϕ)
    x = rand(p.p_x)
    y = A * sin(x + ϕ)
    return x, y
end

# PMAML Sinusoid+Line
struct SinusoidLinearMixture <: AbstractSyntheticDataSet
    p_sinusoid::Bernoulli
    # Sinusoid
    p_A::UnivariateDistribution
    p_ϕ::UnivariateDistribution
    # Line
    p_a::UnivariateDistribution
    p_b::UnivariateDistribution
    # Marginal
    p_x::UnivariateDistribution
    # Error
    p_ϵ::UnivariateDistribution
end
function SinusoidLinearMixture()
    return SinusoidLinearMixture(Bernoulli(),
                                 Uniform(0.1, 5.0),
                                 Uniform(0.0, π),
                                 Uniform(-3.0, 3.0),
                                 Uniform(-3.0, 3.0),
                                 Uniform(-5.0, 5.0),
                                 Normal(0.0, 0.3^2))
end

function Base.rand(p::SinusoidLinearMixture, k::Int)
    X = rand(p.p_x, k)
    if rand(p.p_sinusoid)
        A = rand(p.p_A)
        ϕ = rand(p.p_ϕ)
        Y = @. A * sin(X + ϕ)
    else
        a = rand(p.p_a)
        b = rand(p.p_b)
        Y = @. a * X + b
    end
    return X, Y
end
function Base.rand(p::SinusoidLinearMixture)
    x = rand(p.p_x)
    if rand(p.p_sinusoid)
        A = rand(p.p_A)
        ϕ = rand(p.p_ϕ)
        y = A * sin(X + ϕ)
    else
        a = rand(p.p_a)
        b = rand(p.p_b)
        y = a * X + b
    end
    return x, y
end

# Uncertainty in Multitask Transfer learning, Harmonics
# NOTE: Unclear what sigma_y is from paper
struct Harmonic <: AbstractSyntheticDataSet
    p_A::UnivariateDistribution
    p_ω::UnivariateDistribution
    p_ϕ::UnivariateDistribution
    p_meta_x::UnivariateDistribution
    p_ϵ::UnivariateDistribution
end
Harmonic() = Harmonic(Uniform(0.0, 1.0), Uniform(5.0, 7.0), Uniform(0, 2π), Uniform(-4.0, 4.0), Normal(0.0, 0.1))

function Base.rand(p::Harmonic, k::Int)
    A = rand(p.p_A, 2)
    A₁, A₂ = A[1], A[2]
    ω = rand(p.p_ω)
    ϕ = rand(p.p_ϕ, 2)
    ϕ₁, ϕ₂ = ϕ[1], ϕ[2]
    μ_x = rand(p.p_meta_x)
    p_x = Normal(μ_x, 1.0)
    X = rand(p_x, k)
    ϵ = rand(p.p_ϵ, k)
    Y = @. A₁ * sin(ω * X + ϕ₁) + A₂ * sin(2ω * X + ϕ₂) + ϵ
    return X, Y
end
function Base.rand(p::Harmonic)
    A = rand(p.p_A, 2)
    A₁, A₂ = A[1], A[2]
    ω = rand(p.p_ω)
    ϕ = rand(p.p_ϕ, 2)
    ϕ₁, ϕ₂ = ϕ[1], ϕ[2]
    μ_x = rand(p.p_meta_x)
    p_x = Normal(μ_x, 1.0)
    X = rand(p_x)
    ϵ = rand(p.p_ϵ)
    y = A₁ * sin(ω * x + ϕ₁) + A₂ * sin(2ω * x + ϕ₂) + ϵ
    return x, y
end


text:
module BioSimulator

using Random, DataStructures, Statistics, PoissonRandom
using StaticArrays, SparseArrays, LinearAlgebra
using RecipesBase

## imports
import RecursiveArrayTools: AbstractDiffEqArray

import Combinatorics: multiexponents
import MacroTools: postwalk
import StatsBase: sample
import Base: show, <=
import Base: getindex, setindex!, iterate, firstindex, lastindex, length, isempty, haskey, get, in
import Base: eachindex, push!, empty!, copy
import Base.Order: Ordering, ForwardOrdering, ReverseOrdering, Forward, Reverse

## load order

## interface

#### well-mixed
include(joinpath("interface", "well-mixed", "species.jl"))
include(joinpath("interface", "well-mixed", "reaction.jl"))
include(joinpath("interface", "well-mixed", "network.jl"))

export Species, Reaction, Network
export <=, number_species, number_reactions, species_list, reaction_list
export get_species, get_reaction

#### discrete-spatial
include(joinpath("interface", "spatial-discrete", "reaction.jl"))

# export @def_reactions, @enumerate_with_sclass, @enumerate_with_nclass
export @def_reactions, @enumerate_with_sclass

## data structures
include(joinpath("data-structures", "dep_graph.jl"))

export DGView, DGLazy, DGVector, dependents, rxnrxn_depgraph, spcrxn_depgraph

include(joinpath("data-structures", "priority_queue.jl"))

export PQBinaryHeap, peektop

## state
include(joinpath("state", "site.jl"))
include(joinpath("state", "vonneumann.jl"))
include(joinpath("state", "hexagonal.jl"))
include(joinpath("state", "Lattice.jl"))

export VonNeumann, Hexagonal
export Lattice

const global NBTYPES = [VonNeumann(), Hexagonal()]

## model
include(joinpath("model", "reaction_system.jl"))
include(joinpath("model", "SampleClassEnumeration.jl"))
include(joinpath("model", "ips.jl"))

## algorithms
include(joinpath("algorithms", "ssa_utils.jl"))
include(joinpath("algorithms", "abstract_algorithm.jl"))
include(joinpath("algorithms", "direct.jl"))
include(joinpath("algorithms", "enhanced_direct.jl"))
include(joinpath("algorithms", "sorting_direct.jl"))
include(joinpath("algorithms", "firstreaction.jl"))
include(joinpath("algorithms", "nextreaction.jl"))
include(joinpath("algorithms", "rejection.jl"))

include(joinpath("algorithms", "tauleap_utils.jl"))
include(joinpath("algorithms", "poisson_tauleap.jl"))
include(joinpath("algorithms", "stepanticipation.jl"))

export HasRates, HasSums

## simulators
include(joinpath("simulators", "abstract_simulator.jl"))
include(joinpath("simulators", "exact.jl"))
include(joinpath("simulators", "tauleap.jl"))
include(joinpath("simulators", "tauleap_hybrid.jl"))
include(joinpath("simulators", "build_simulator.jl"))

export Direct, EnhancedDirect, SortingDirect
export FirstReaction, NextReaction, RejectionSSA
export TauLeapingDG2001, TauLeapingDGLP2003, StepAnticipation
export HybridSAL

export cumulative_intensity, jump_rates, next_jump_index, next_jump_time

## simulation + output
include("simulate.jl")
include(joinpath("output", "SamplePath.jl"))
include(joinpath("output", "Configuration.jl"))
include(joinpath("output", "tabletraits.jl"))

export simulate, parse_model, SamplePath, Ensemble, Configuration
export save_state, save_rates

end # module

ext: jl


text:
using OffsetArrays, LinearAlgebra, LoopVectorization, Test

let nt = min(Threads.nthreads(), (Sys.CPU_THREADS)::Int)
  if (LoopVectorization.num_cores() < 4) && (nt ≥ 4)
    @eval LoopVectorization.num_cores() = LoopVectorization.StaticInt{$nt}()
  end
end

function mydotavx(a, b)
  s = zero(eltype(a))
  @tturbo for i ∈ eachindex(a, b)
    s += a[i] * b[i]
  end
  s
end
function AmulB!(C, A, B)
  @tturbo for n in indices((B, C), 2), m in indices((A, C), 1)
    Cmn = zero(eltype(C))
    for k in indices((A, B), (2, 1))
      Cmn += A[m, k] * B[k, n]
    end
    C[m, n] = Cmn
  end
  C
end
function dot3(
  x::AbstractVector{Complex{T}},
  A::AbstractMatrix{Complex{T}},
  y::AbstractVector{Complex{T}},
) where {T}
  xr = reinterpret(reshape, T, x)
  yr = reinterpret(reshape, T, y)
  Ar = reinterpret(reshape, T, A)
  sre = zero(T)
  sim = zero(T)
  @tturbo for n in axes(Ar, 3)
    tre = zero(T)
    tim = zero(T)
    for m in axes(Ar, 2)
      tre += xr[1, m] * Ar[1, m, n] + xr[2, m] * Ar[2, m, n]
      tim += xr[1, m] * Ar[2, m, n] - xr[2, m] * Ar[1, m, n]
    end
    sre += tre * yr[1, n] - tim * yr[2, n]
    sim += tre * yr[2, n] + tim * yr[1, n]
  end
  sre + im * sim
end
function conv!(out, A, kern)
  @tturbo for n ∈ axes(out, 2), m ∈ axes(out, 1)
    tmp = zero(eltype(out))
    for nᵢ ∈ -1:1, mᵢ ∈ -1:1
      tmp += A[mᵢ+m, nᵢ+n] * kern[mᵢ, nᵢ]
    end
    out[m, n] = tmp
  end
  out
end
function conv_baseline!(out, A, kern)
  for n ∈ axes(out, 2), m ∈ axes(out, 1)
    tmp = zero(eltype(out))
    for nᵢ ∈ -1:1, mᵢ ∈ -1:1
      tmp += A[mᵢ+m, nᵢ+n] * kern[mᵢ, nᵢ]
    end
    out[m, n] = tmp
  end
  out
end

struct DenseConvDims{N,K,C_in,C_out} end

function kernaxes(::DenseConvDims{2,K,C_in,C_out}) where {K,C_in,C_out}
  K₁ = LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(K[1])
  K₂ = LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(K[2])
  Cᵢₙ = LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(C_in)
  Cₒᵤₜ = LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(C_out)
  (K₁, K₂, Cᵢₙ, Cₒᵤₜ)
end

function convlayer!(
  out::AbstractArray{<:Any,4},
  img,
  kern,
  dcd::DenseConvDims{2,<:Any,<:Any,<:Any},
)
  (K₁, K₂, Cᵢₙ, Cₒᵤₜ) = kernaxes(dcd)
  @tturbo for j₁ ∈ axes(out, 1), j₂ ∈ axes(out, 2), d ∈ axes(out, 4), o ∈ Cₒᵤₜ
    s = zero(eltype(out))
    for k₁ ∈ K₁, k₂ ∈ K₂, i ∈ Cᵢₙ
      s += img[j₁+k₁-1, j₂+k₂-1, i, d] * kern[k₁, k₂, i, o]
    end
    out[j₁, j₂, o, d] = s
  end
  out
end
function convlayer_direct!(
  out::AbstractArray{<:Any,4},
  img,
  kern,
  dcd::DenseConvDims{2,<:Any,<:Any,<:Any},
)
  (K₁, K₂, Cᵢₙ, Cₒᵤₜ) = kernaxes(dcd)
  @inbounds @fastmath for j₁ ∈ axes(out, 1), j₂ ∈ axes(out, 2), d ∈ axes(out, 4), o ∈ Cₒᵤₜ
    s = zero(eltype(out))
    for k₁ ∈ K₁, k₂ ∈ K₂, i ∈ Cᵢₙ
      s += img[j₁+k₁-1, j₂+k₂-1, i, d] * kern[k₁, k₂, i, o]
    end
    out[j₁, j₂, o, d] = s
  end
  out
end

function thread_cartesian_indices!(dst, src)
  @turbo thread = true for i in eachindex(dst, src)
    dst[i] = src[i]
  end
end

@testset "Threading" begin
  @show @__LINE__
  dcd = DenseConvDims{2,(5, 5),3,6}()
  kern4 = rand(Float32, 5, 5, 3, 6)
  @time for M ∈ 17:50:267
    img = rand(Float32, M, M, 3, 100)
    outimage1 = Array{Float32}(
      undef,
      size(img, 1) + 1 - size(kern4, 1),
      size(img, 2) + 1 - size(kern4, 2),
      size(kern4, 4),
      size(img, 4),
    )
    outimage2 = similar(outimage1)

    convlayer!(outimage1, img, kern4, dcd)
    convlayer_direct!(outimage2, img, kern4, dcd)
    @test outimage1 ≈ outimage2
  end

  @time for M ∈ 17:399
    # @show M
    K = M
    N = M
    A = rand(M, K)
    B = rand(K, N)
    b = rand(K)
    @test dot(A, B) ≈ mydotavx(A, B)

    C1 = A * B
    @test AmulB!(similar(C1), A, B) ≈ C1
    c1 = A * b
    @test AmulB!(similar(c1), A, b) ≈ c1

    if VERSION ≥ v"1.6"
      x = randn(Complex{Float64}, 3M - 1)
      W = randn(Complex{Float64}, 3M - 1, 3M + 1)
      y = randn(Complex{Float64}, 3M + 1)
      @test dot(x, W, y) ≈ dot3(x, W, y)
    end

    kern = OffsetArray(randn(3, 3), -2, -2)
    out1 = OffsetArray(randn(size(A) .- 2), 1, 1)
    out2 = similar(out1)
    @test conv!(out1, A, kern) ≈ conv_baseline!(out2, A, kern)
  end

  src_big = rand(4, 10)
  dst_big = similar(src_big)
  src = view(src_big, 2:3, :)
  dst = view(dst_big, 2:3, :)
  thread_cartesian_indices!(dst, src)
  dst ≈ src

end


text:
# set discrete values of mu for input observations
const disc_mu = [0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.85, 0.9, 0.95, 1.0]
const mu_symb = [:mu02, :mu03, :mu04, :mu05, :mu06, :mu07, :mu08, :mu085, :mu09, :mu095, :mu10]
const disc_ax = [:n, :e, :s, :w, :c]

make_grid(N::Integer) = range(-1.0, 1.0, length=N)
make_grid(;N::Integer=256) = make_grid(N)

function calc_r2(x::T,y::T) where T<:AF
    return x*x + y*y
end

function calc_r2(t::Tuple{T,T}) where T<:AF
    return calc_r2(t[1], t[2])
end

# Get mu on disk from provided (x,y) position:
# TODO: relies on small angle approx?
function calc_mu(t::Tuple{T,T}) where T<:AF
    r2 = calc_r2(t)
    return sqrt((r2 < one(T)) * (one(T) - r2))
end

function calc_mu(x::T, y::T) where T<:AF
    return calc_mu((x,y))
end

# Calculate mu for each position on a grid
function mu_map(grid::AA{T,1}) where T<:AF
    return calc_mu.((x,y) for x in grid, y in grid)
end

function mu_map(N::Integer)
    return mu_map(make_grid(N))
end

function calc_norm_term(x::T, y::T, N::Integer, u1::T, u2::T) where T<:AF
    return quad_limb_darkening(x,y,u1,u2) * π / (2.0 * N^2)
end

function calc_norm_term(x::T, y::T, disk::DiskParams) where T<:AF
    return calc_norm_term(x, y, disk.N, disk.u1, disk.u2)
end

# Find the nearest axis to a given point on a grid
function find_nearest_ax(x::T, y::T) where T<:AF
    if (x^2 + y^2) > one(T)
        return :off
    elseif ((y == zero(T)) & (x == zero(T))) # center
        return :c
    elseif y >= abs(x) # north
        return :n
    elseif y <= -abs(x) # south
        return :s
    elseif x <= -abs(y) # east
        return :e
    elseif x >= abs(y) # west
        return :w
    end
end

function find_nearest_mu(x::T, y::T) where T<:AF
    return searchsortednearest(disc_mu, calc_mu(x, y))
end

function assemble_dict_key(mu_ind::Int, ax::Symbol; mu_symb::AA{Symbol,1}=mu_symb)
    if mu_ind == 11 # e.g., if mu == 1.0
        return (:c, :mu10)
    else
        return (ax, mu_symb[mu_ind])
    end
end

function get_key_for_pos(x::T, y::T) where T<:AF
    # find nearest mu index
    mu_ind = find_nearest_mu(x, y)

    # get nearest axis and mu
    near_mu = disc_mu[mu_ind]
    near_ax = find_nearest_ax(x, y)

    # assemble dictionary key & extract values
    key = assemble_dict_key(mu_ind, near_ax)
    return key
end

function mu_to_xy(mu::T, ax::Symbol) where T<:AF
    @assert zero(T) <= mu <= one(T)
    @assert ax in [:n, :e, :s, :w, :c]

    # get radial distance from disk center
    r = sqrt(one(T) - mu^2)

    if ax == :n
        return (0.0, r)
    elseif ax == :e
        return (-r, 0.0)
    elseif ax == :w
        return (r, 0.0)
    elseif ax == :s
        return (0.0, -r)
    else
        return (0.0, 0.0)
    end
end

max_stars_repo_name: palumbom/GRASS

max_issues_repo_name: palumbom/GRASS

max_forks_repo_name: palumbom/GRASS


text:
using FastaIO;
using DataFrames;

"""
parse_input(seq_fp1::String, seq_fp2::String, sm_fp::String)
Read in FASTA files from SEQ_FP1 and SEQ_FP2 and the substitution matrix
stored in SM_FP.

"""
function parse_input(seq_fp1, seq_fp2, sm_fp)

    seq1 = readfasta(seq_fp1)[1][2];
    seq2 = readfasta(seq_fp2)[1][2];

    sm = parse_score_matrix(sm_fp);

    seq1, seq2, sm
end

"""
parse_score_matrix(sm_fp::String)
Read in subsitution matrix from SM_FP and return a DataFrame.
"""
function parse_score_matrix(sm_fp)

    f = open(sm_fp);
    lines = readlines(f);
    i = 1
    while lines[i][1] == '#'
        i += 1
    end
    #aas = strip.(split(lines[i], "  "));
    aas = strip.(split(lines[i]))
    aas = [convert(String, aa) for aa in aas];
    i += 1

    N = length(aas);

    #sm = DataFrame(zeros(N));
    sm = DataFrame(Float64, 0, N);
    nframe = DataFrame(AA = aas);

    while i <= length(lines)

        sub = strip.(split(lines[i]));
        sub = parse.(Float64, sub);
        push!(sm, sub);
        i += 1
    end

    close(f)
    #deleterows!(sm, 1);
    names!(sm, [Symbol(aa) for aa in aas]);
    sm = hcat(nframe, sm);

    sm

end
"""
parse_list(seqs_fp::String)
Read in a set of sequences to be aligned, store in an Array.
"""
function parse_list(seqs_fp)

    seqs = []

    f = open(seqs_fp)
    lines = readlines(f)

    for l in lines
        seq_pair = strip.(split(l));
        seq_pair = convert.(String, seq_pair);
        seq1 = readfasta(seq_pair[1])[1][2]
        seq2 = readfasta(seq_pair[2])[1][2];
        push!(seqs, [seq1, seq2])
    end

    close(f)

    seqs

end
"""
parse_seq_list(seqs_fp::String, sm_fp::String)
Read in a set of sequences to be aligned, store in an Array; read in the
substitution score matrix stored at SM_FP as well. Parsing function for some analyses.
"""
function parse_seq_list(seqs_fp, sm_fp)

    seqs = parse_list(seqs_fp);
    sm = parse_score_matrix(sm_fp);

    seqs, sm

end

"""
parse_for_benchmarking(pos_seqs_fp, neg_seqs_fp, sm_fp)
Read in set of positive sequence pairs POS_SEQS_FP and negative sequence pairs
NEG_SEQS_FP for ROC benchmarking. Also read in the score matrix stored in SM_FP.
"""
function parse_for_benchmarking(pos_seqs_fp, neg_seqs_fp, sm_fp)

    pos_seqs = parse_list(pos_seqs_fp)
    neg_seqs = parse_list(neg_seqs_fp)
    sm = parse_score_matrix(sm_fp)

    pos_seqs, neg_seqs, sm

end


text:
using LaTeXStrings
using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity

import ParameterEstimocean.Parameters: closure_with_parameters

# catke_parameter_guide = Dict(:Cᴰ => (name = "Dissipation parameter (TKE equation)", latex = L"C^D", default = 2.9079, bounds = (0.0, 0.5)),
#                         :Cᴸᵇ => (name = "Mixing length parameter", latex = L"C^{\ell}_b", default = 1.1612, bounds = (0.0, 10.0)),
#                         :Cᵟu => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}u", default = 0.5, bounds = (0.0, 3.0)),
#                         :Cᵟc => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}c", default = 0.5, bounds = (0.0, 3.0)),
#                         :Cᵟe => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}e", default = 0.5, bounds = (0.0, 3.0)),
#                         :Cᵂu★ => (name = "Mixing length parameter", latex = L"C^W_{u\star}", default = 3.6188, bounds = (0.0, 1.0)),
#                         :CᵂwΔ => (name = "Mixing length parameter", latex = L"C^Ww\Delta", default = 1.3052, bounds = (0.0, 30.0)),
#                         :CᴷRiʷ => (name = "Stability function parameter", latex = L"C^KRi^w", default = 0.7213, bounds = (0.0, 20.0)),
#                         :CᴷRiᶜ => (name = "Stability function parameter", latex = L"C^KRi^c", default = 0.7588, bounds = (-1.5, 3.0)),
#                         :Cᴷu⁻ => (name = "Velocity diffusivity LB", latex = L"C^Ku^-", default = 0.1513, bounds = (0.0, 0.5)),
#                         :Cᴷuʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Ku^r", default = 3.8493, bounds = (0.0, 3.0)),
#                         :Cᴷc⁻ => (name = "Tracer diffusivity LB", latex = L"C^Kc^-", default = 0.3977, bounds = (0.0, 0.8)),
#                         :Cᴷcʳ => (name = "Tracer diffusivity (UB-LB)/LB", latex = L"C^Kc^r", default = 3.4601, bounds = (0.0, 30.0)),
#                         :Cᴷe⁻ => (name = "TKE diffusivity LB", latex = L"C^Ke^-", default = 0.1334, bounds = (0.0, 0.5)),
#                         :Cᴷeʳ => (name = "TKE diffusivity (UB-LB)/LB", latex = L"C^Ke^r", default = 8.1806, bounds = (0.0, 60.0)),
#                         :Cᴬu => (name = "Convective adjustment velocity parameter", latex = L"C^A_U", default = 0.0057, bounds = (0.0, 0.2)),
#                         :Cᴬc => (name = "Convective adjustment tracer parameter", latex = L"C^A_C", default = 0.6706, bounds = (0.0, 2.0)),
#                         :Cᴬe => (name = "Convective adjustment TKE parameter", latex = L"C^A_E", default = 0.2717, bounds = (0.0, 2.0)),
#                         :Cᵇu => (name = "Convective adjustment TKE parameter", latex = L"C^b_U", default = 0.596, bounds = (0.0, 1.5)),
#                         :Cᵇc => (name = "Convective adjustment TKE parameter", latex = L"C^b_C", default = 0.0723, bounds = (0.0, 0.3)),
#                         :Cᵇe => (name = "Convective adjustment TKE parameter", latex = L"C^b_E", default = 0.637, bounds = (0.0, 0.6)),
#                         :Cˢu => (name = "Convective adjustment TKE parameter", latex = L"C^s_U", default = 0.628, bounds = (0.0, 4.0)),
#                         :Cˢc => (name = "Convective adjustment TKE parameter", latex = L"C^s_C", default = 0.426, bounds = (0.0, 4.0)),
#                         :Cˢe => (name = "Convective adjustment TKE parameter", latex = L"C^s_E", default = 0.711, bounds = (0.0, 2.0)),
# )
catke_parameter_guide = Dict(:Cᴰ => (name = "Dissipation parameter (TKE equation)", latex = L"C^D", default = 2.9079, bounds = (0.0, 1.5)),
                        :Cᴸᵇ => (name = "Mixing length parameter", latex = L"C^{\ell}_b", default = 1.1612, bounds = (0.0, 10.0)),
                        :Cᵟu => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}u", default = 0.5, bounds = (0.0, 3.0)),
                        :Cᵟc => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}c", default = 0.5, bounds = (0.0, 3.0)),
                        :Cᵟe => (name = "Ratio of mixing length to grid spacing", latex = L"C^{\delta}e", default = 0.5, bounds = (0.0, 3.0)),
                        :Cᵂu★ => (name = "Mixing length parameter", latex = L"C^W_{u\star}", default = 3.6188, bounds = (0.0, 2.5)),
                        :CᵂwΔ => (name = "Mixing length parameter", latex = L"C^Ww\Delta", default = 1.3052, bounds = (0.0, 70.0)),
                        :CᴷRiʷ => (name = "Stability function parameter", latex = L"C^KRi^w", default = 0.7213, bounds = (0.0, 50.0)),
                        :CᴷRiᶜ => (name = "Stability function parameter", latex = L"C^KRi^c", default = 0.7588, bounds = (-1.5, 3.0)),
                        :Cᴷu⁻ => (name = "Velocity diffusivity LB", latex = L"C^Ku^-", default = 0.1513, bounds = (0.0, 0.5)),
                        :Cᴷuʳ => (name = "Velocity diffusivity (UB-LB)/LB", latex = L"C^Ku^r", default = 3.8493, bounds = (0.0, 6.0)),
                        :Cᴷc⁻ => (name = "Tracer diffusivity LB", latex = L"C^Kc^-", default = 0.3977, bounds = (0.0, 1.2)),
                        :Cᴷcʳ => (name = "Tracer diffusivity (UB-LB)/LB", latex = L"C^Kc^r", default = 3.4601, bounds = (0.0, 40.0)),
                        :Cᴷe⁻ => (name = "TKE diffusivity LB", latex = L"C^Ke^-", default = 0.1334, bounds = (0.0, 0.8)),
                        :Cᴷeʳ => (name = "TKE diffusivity (UB-LB)/LB", latex = L"C^Ke^r", default = 8.1806, bounds = (0.0, 100.0)),
                        :Cᴬu => (name = "Convective adjustment velocity parameter", latex = L"C^A_U", default = 0.0057, bounds = (0.0, 0.2)),
                        :Cᴬc => (name = "Convective adjustment tracer parameter", latex = L"C^A_C", default = 0.6706, bounds = (0.0, 2.0)),
                        :Cᴬe => (name = "Convective adjustment TKE parameter", latex = L"C^A_E", default = 0.2717, bounds = (0.0, 2.0)),
                        :Cᵇu => (name = "Convective adjustment TKE parameter", latex = L"C^b_U", default = 0.596, bounds = (0.0, 5.0)),
                        :Cᵇc => (name = "Convective adjustment TKE parameter", latex = L"C^b_C", default = 0.0723, bounds = (0.0, 0.6)),
                        :Cᵇe => (name = "Convective adjustment TKE parameter", latex = L"C^b_E", default = 0.637, bounds = (0.0, 1.0)),
                        :Cˢu => (name = "Convective adjustment TKE parameter", latex = L"C^s_U", default = 0.628, bounds = (0.0, 10.0)),
                        :Cˢc => (name = "Convective adjustment TKE parameter", latex = L"C^s_C", default = 0.426, bounds = (0.0, 7.0)),
                        :Cˢe => (name = "Convective adjustment TKE parameter", latex = L"C^s_E", default = 0.711, bounds = (0.0, 2.5)),
)

ri_based_parameter_guide = Dict(:ν₀ => (name = "base viscosity", latex = L"\nu_0", default = 0.01, bounds = (0.0, 1.0)),
                        :κ₀ => (name = "base diffusivity", latex = L"\kappa_0", default = 0.1, bounds = (0.0, 10.0)),
                        :Ri₀ν => (name = "Ri viscosity", latex = L"Ri_0\nu", default = -0.5, bounds = (-20.0, 20.0)),
                        :Ri₀κ => (name = "Ri diffusivity", latex = L"Ri_0\kappa", default = -0.5, bounds = (-20.0, 20.0)),
                        :Riᵟν => (name = "Ri delta viscosity", latex = L"Ri^{\delta}\nu", default = 1.0, bounds = (0.0, 10.0)),
                        :Riᵟκ => (name = "Ri delta diffusivity", latex = L"Ri^{\delta}\kappa", default = 1.0, bounds = (0.0, 3.0)),
)

"""
    ParameterSet{C, N, S} where C

Parameter set containing the names `names` of parameters, and a
NamedTuple `settings` mapping names of "background" parameters
to their fixed values to be maintained throughout the calibration.
"""
struct ParameterSet{C, N, S}
       names :: N
    settings :: S

    function ParameterSet{C}(names::N = nothing, settings::S = nothing) where {C, N, S}
        return new{C, N, S}(names, settings)
    end
end

parameter_guide(::ParameterSet{<:CATKEVerticalDiffusivity}) = catke_parameter_guide
parameter_guide(::ParameterSet{<:RiBasedVerticalDiffusivity}) = ri_based_parameter_guide

bounds(name, parameter_set) = parameter_guide(parameter_set)[name].bounds
default(name, parameter_set) = parameter_guide(parameter_set)[name].default

function named_tuple_map(names, f)
    names = Tuple(names)
    return NamedTuple{names}(f.(names))
end

"""
    ParameterSet(names::Set; nullify = Set())

Construct a `ParameterSet` containing all of the information necessary
to build a closure with the specified default parameters and settings,
given a Set `names` of the parameter names to be tuned, and a Set `nullify`
of parameters to be set to zero.
"""
function ParameterSet{C}(names::Set; nullify = Set()) where C
    ref_set = ParameterSet{C}()
    zero_set = named_tuple_map(nullify, name -> 0.0)
    bkgd_set = named_tuple_map(keys(parameter_guide(ref_set)), name -> default(name, ref_set))
    settings = merge(bkgd_set, zero_set) # order matters: `zero_set` overrides `bkgd_set`
    return ParameterSet{C}(Tuple(names), settings)
end

names(ps::ParameterSet) = ps.names

###
### Define some convenient parameter sets
###

# CATKEVerticalDiffusivity

required_params = Set([:Cᵟu, :Cᵟc, :Cᵟe, :Cᴰ, :Cᴸᵇ, :Cᵂu★, :CᵂwΔ, :Cᴷu⁻, :Cᴷc⁻, :Cᴷe⁻])
ri_depen_params = Set([:CᴷRiʷ, :CᴷRiᶜ, :Cᴷuʳ, :Cᴷcʳ, :Cᴷeʳ])
conv_adj_params = Set([:Cᴬu, :Cᴬc, :Cᴬe])

CATKEParametersRiDependent = ParameterSet{CATKEVerticalDiffusivity}(union(required_params, ri_depen_params); nullify = conv_adj_params)
CATKEParametersRiIndependent = ParameterSet{CATKEVerticalDiffusivity}(union(required_params); nullify = union(conv_adj_params, ri_depen_params))
CATKEParametersRiDependentConvectiveAdjustment = ParameterSet{CATKEVerticalDiffusivity}(union(required_params, conv_adj_params, ri_depen_params))
CATKEParametersRiIndependentConvectiveAdjustment = ParameterSet{CATKEVerticalDiffusivity}(union(required_params, conv_adj_params); nullify = ri_depen_params)

# RiBasedVerticalDiffusivity
ri_based_params = Set([:ν₀, :κ₀, :Ri₀ν, :Ri₀κ, :Riᵟν, :Riᵟκ])

RiBasedParameterSet = ParameterSet{RiBasedVerticalDiffusivity}(ri_based_params)


text:
using HistogramThresholding, TestImages, ColorTypes, ImageContrastAdjustment
using StaticArrays
using Test, Aqua, Documenter

@testset "Project meta quality checks" begin
    # Not checking compat section for test-only dependencies
    Aqua.test_all(HistogramThresholding; project_extras=true, deps_compat=true, stale_deps=true, project_toml_formatting=true)

    DocMeta.setdocmeta!(HistogramThresholding, :DocTestSetup, :(using HistogramThresholding); recursive=true)
    doctest(HistogramThresholding, manual = false)
end

function gaussian_mixture(z;P₁,μ₁,σ₁,P₂,μ₂,σ₂)
    (P₁ / (σ₁*√(2*π))) * exp( -(z-μ₁)^2 / (2*σ₁^2) )  +  (P₂ / (σ₂*√(2*π))) * exp( -(z-μ₂)^2 / (2*σ₂^2) )
end

function surge(t;A,k)
    A*t*exp(-k*t)
end

@testset "HistogramThresholding.jl" begin
    include("otsu.jl")
    include("minimum.jl")
    include("intermodes.jl")
    include("unimodal.jl")
    include("moments.jl")
    include("minimum_error.jl")
    include("balancedthreshold.jl")
    include("yen.jl")
    include("entropy_thresholding.jl")
end


text:
# module Construction

# using LightGraphs, MetaGraphs
# using ..Arrays

export
    MatrixWrapper,
    IndicatorGrid,
    VtxGrid

abstract type MatrixWrapper <: AbstractMatrix{Int} end
for op in [:getindex,:setindex!,:get,:get!,:size,:length,:iterate]
    @eval Base.$op(m::MatrixWrapper,args...) = $op(m.mtx,args...)
end
Base.get(m::MatrixWrapper,idxs::NTuple{N,Int},args...) where {N} = get(m.mtx,idxs,args...)

"""
    IndicatorGrid

A matrix such that `M[i,j] == 0` indicates "free" and `M[i,j] == 1` indicates
"obstacle".
"""
struct IndicatorGrid{M<:AbstractMatrix} <: MatrixWrapper
    mtx::M
end
function Base.getindex(m::M,ind1::I,ind2::J) where {M<:IndicatorGrid,I<:AbstractRange{Int},J<:AbstractRange{Int}}
    M(getindex(m.mtx,ind1,ind2))
end

"""
    VtxGrid

A matrix such that `M[i,j] == v`, where v is the index of the
vertex whose coordinates are (i,j)
"""
struct VtxGrid{M<:AbstractMatrix} <: MatrixWrapper
    mtx::M
end

for t in [:IndicatorGrid,:VtxGrid]
    @eval $t(m::$t) = $t(m.mtx)
end

export
    initialize_vtx_grid_from_indicator_grid,
    initialize_regular_vtx_grid,
    initialize_dense_vtx_grid,
    construct_vtx_grid,
    initialize_grid_graph_from_vtx_grid,
    initialize_regular_grid_graph,
    initialize_grid_graph_with_obstacles

"""
    initialize_vtx_grid_from_indicator_grid()

Args:
    * grid: an indicator grid, where grid[i,j] == 0 denotes that cell (i,j) is
    free, and grid[i,j] != 0 denotes that the cell is impassable.

Returns:
    * A matrix whose non-zero elements represent the ids of vertices in a graph
    embedded in the 2D grid.
"""
function initialize_vtx_grid_from_indicator_grid(grid)
    K = VtxGrid(zeros(Int,size(grid)))
    idx = 1
    for i in 1:size(K,1)
        for j in 1:size(K,2)
            if grid[i,j] == 0
                K[i,j] = idx
                idx += 1
            end
        end
    end
    K
end
VtxGrid(m::IndicatorGrid) = initialize_vtx_grid_from_indicator_grid(m)

construct_indicator_grid_from_vtx_grid(vtx_grid) = IndicatorGrid(Int.(vtx_grid .== 0))
IndicatorGrid(m::VtxGrid) = construct_indicator_grid_from_vtx_grid(m)

export convolve_with_occupancy_kernel

"""
    convolve_with_occupancy_kernel(indicator_grid::IndicatorGrid,kernel)

Construct a new indicator grid corresponding to convolution of the original grid
with a kernel. The second argument may be a Matrix{Int} or a Tuple{Int,Int}
indicating the size for a ones() kernel.
"""
function convolve_with_occupancy_kernel(indicator_grid::IndicatorGrid,kernel::Matrix{Int})
    filtered_grid = IndicatorGrid(
        Int.(imfilter(indicator_grid, centered(kernel),Inner()) .> 0)
        )
end
function convolve_with_occupancy_kernel(indicator_grid::IndicatorGrid,shape::Tuple{Int,Int})
    convolve_with_occupancy_kernel(indicator_grid,ones(Int,shape))
end

export
    vtx_list_from_vtx_grid

function vtx_list_from_vtx_grid(m::VtxGrid)
    vtx_list = Vector{Tuple{Int,Int}}()
    for i in 1:size(m,1)
        for j in 1:size(m,2)
            if m[i,j] > 0
                push!(vtx_list,(i,j))
                @assert length(vtx_list) == m[i,j]
            end
        end
    end
    vtx_list
end

"""
    Returns a grid that represents a 2D environment with regularly spaced
    rectangular obstacles.

    Indicator grid values:
        0 == FREE
        1 == OCCUPIED
"""
function initialize_regular_indicator_grid(;
    n_obstacles_x=2,
    n_obstacles_y=2,
    obs_width = [2;2],
    obs_offset = [1;1],
    env_pad = [1;1]
    )
    # generate occupancy grid representing the environment
    o = ones(Int,obs_width[1],obs_width[2]) # obstacle region
    op = pad_matrix(o,(obs_offset[1],obs_offset[2]),0) # padded obstacles region
    A = repeat(op,n_obstacles_x,n_obstacles_y)
    Ap = pad_matrix(A,(env_pad[1],env_pad[2]),0) # padded occupancy grid
    IndicatorGrid(Ap)
end
function initialize_regular_vtx_grid(;kwargs...)
    grid = initialize_regular_indicator_grid(;kwargs...)
    initialize_vtx_grid_from_indicator_grid(grid)
end

"""
    `initialize_dense_vtx_grid`
"""
initialize_dense_vtx_grid(x_dim::Int,y_dim::Int) = VtxGrid(collect(reshape(collect(1:x_dim*y_dim),y_dim,x_dim)'))

"""
    `construct_vtx_grid`

Returns a matrix `M` such that `MV[i,j] == v`, where v is the index of the
vertex whose coordinates are (i,j)

Arguments:
* vtxs : a list of integer coordinates
* dims : the dimensions of the grid
"""
function construct_vtx_grid(vtxs::V,dims) where {V<:Vector}
    vtx_grid = VtxGrid(zeros(Int,dims))
    for (k,vtx) in enumerate(vtxs)
        vtx_grid[vtx...] = k
    end
    vtx_grid
end

"""
    Returns a grid graph that represents a 2D environment with regularly spaced
    rectangular obstacles
"""
function initialize_grid_graph_from_vtx_grid(K::M) where {M<:Union{VtxGrid,Matrix{Int}}}
    G = MetaGraph()
    for i in 1:size(K,1)
        for j in 1:size(K,2)
            if K[i,j] != 0
                add_vertex!(G)
                add_edge!(G,nv(G),nv(G))
            end
        end
    end

    for i in 1:size(K,1)
        for j in 1:size(K,2)
            if K[i,j] != 0
                if j < size(K,2)
                    add_edge!(G,K[i,j],K[i,j+1],Dict(:direction=>:NORTH))
                end
                if j > 1
                    add_edge!(G,K[i,j],K[i,j-1],Dict(:direction=>:SOUTH))
                end
                if i < size(K,1)
                    add_edge!(G,K[i,j],K[i+1,j],Dict(:direction=>:EAST))
                end
                if i > 1
                    add_edge!(G,K[i,j],K[i-1,j],Dict(:direction=>:WEST))
                end
            end
        end
    end
    G
end

"""
    Returns a grid graph that represents a 2D environment with regularly spaced
    rectangular obstacles
"""
function initialize_regular_grid_graph(;
    env_offset = [1.0,1.0],
    env_scale = 2.0, # this is essentially the robot diameter
    kwargs...)

    K = initialize_regular_vtx_grid(;kwargs...)

    graph = initialize_grid_graph_from_vtx_grid(K)
    for i in 1:size(K,1)
        for j in 1:size(K,2)
            if K[i,j] != 0
                set_props!(graph, K[i,j],
                    Dict(:x=>env_offset[1] + env_scale*(i-1),
                    :y=>env_offset[2] + env_scale*(j-1))
                    )
            end
        end
    end
    return graph
end

"""
    Returns a grid graph that represents a 2D environment with obstacles placed
    over specific vertices
"""
function initialize_grid_graph_with_obstacles(
    dims,obstacles::Vector{Tuple{Int,Int}}=Vector{Tuple{Int,Int}}())
    indicator_grid = IndicatorGrid(zeros(Int,dims))
    for obs in obstacles
        indicator_grid[obs...] = 0
    end
    vtx_grid = VtxGrid(indicator_grid)
    G = initialize_grid_graph_from_vtx_grid(vtx_grid)
    return G, vtx_grid
end

# end


text:
"""
A [MIPLIB 2017](http://miplib.zib.de/index.html) instance stored in a
`MathOptInterface` backing model.

    MIPLIBInstance(instance_name::String)

Pulls down the MIPLIB instance named `instance_name` from the official
MIPLIB website using `curl`.

Some pieces of the [functional interface](@ref functional_interface) are
implemented:

    prep_instance!(instance::MIPLIBInstance)

!!! note
    You must call `prep_instance!` on a `MIPLIBInstance` to populate the
    internal MathOptInterface model backing it.
"""
mutable struct MIPLIBInstance <: AbstractMOIBackedInstance
    instance_name::String
    tmp_file::String
    model::MOI.FileFormats.MPS.Model

    function MIPLIBInstance(instance_name::String)
        return new(instance_name, tempname())
    end
end

function get_moi_model(inst::MIPLIBInstance)
    @assert isdefined(inst, :model)
    return inst.model
end

const MIPLIB_URL_PREFIX = "http://miplib.zib.de/WebData/instances/"

function prep_instance!(instance::MIPLIBInstance)
    miplib_url = string(MIPLIB_URL_PREFIX, instance.instance_name, ".mps.gz")
    run(`curl $(miplib_url) -o $(instance.tmp_file)`)
    io = GZip.open(instance.tmp_file)
    model_mps = MOI.FileFormats.MPS.Model()
    read!(io, model_mps)
    instance.model = model_mps
    return nothing
end


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIT = MOI.Test

const BRIDGED = MOIB.full_bridge_optimizer(SCIP.Optimizer(display_verblevel=0), Float64)
const CONFIG = MOIT.TestConfig(atol=1e-5, rtol=1e-5, duals=false, infeas_certificates=false)
const CONFIG3 = MOIT.TestConfig(atol=1e-3, rtol=1e-2, duals=false, infeas_certificates=false)

@testset "MOI Continuous Linear" begin
    excluded = [
        "linear1",  # needs MOI.delete (of variables in constraints)
        "linear5",  # needs MOI.delete (of variables in constraints)
        "linear11", # broken in SCIP (#2556)
        "linear13", # TODO: support MOI.FEASIBILITY_SENSE
        "linear14", # needs MOI.delete (of variables in constraints)
    ]
    MOIT.contlineartest(BRIDGED, CONFIG, excluded)
end

@testset "MOI Continuous Conic" begin
    MOIT.lintest(BRIDGED, CONFIG)

    # SOC tests fail because of lower bound requirement of RHS var.
    # MOIT.soctest(BRIDGED, CONFIG)
    # MOIT.rsoctest(BRIDGED, CONFIG)

    # other cones not supported
    # MOIT.geomeantest(BRIDGED, CONFIG)
    # MOIT.exptest(BRIDGED, CONFIG)
    # MOIT.sdptest(BRIDGED, CONFIG)
    # MOIT.logdettest(BRIDGED, CONFIG)
    # MOIT.rootdettest(BRIDGED, CONFIG)
end

@testset "MOI Quadratic Constraint" begin
    # needs objective bridge (MOI/#529)
    # MOIT.qptest(BRIDGED, CONFIG)

    MOIT.qcptest(BRIDGED, CONFIG)
    MOIT.socptest(BRIDGED, CONFIG3)
end

@testset "MOI Integer Linear" begin
    excluded = String[]
    MOIT.intlineartest(BRIDGED, CONFIG, excluded)
    MOIT.indicator3_test(BRIDGED, CONFIG)
end

@testset "MOI Integer Conic" begin
    # SOC tests fail because of lower bound requirement of RHS var.
    # MOIT.intconictest(BRIDGED, CONFIG)
end

@testset "MOI NLP" begin
    # None of tests provide expression graphs in the evaluator.
    # MOIT.nlptest(BRIDGED, CONFIG)
end

@testset "MOI Unit tests" begin
    # TODO: most tests need get-variable-by-name etc.
    # MOIT.unittest(BRIDGED, CONFIG)
end


text:
################################################################################
#                          TEST: DCM <=> Euler Angles
################################################################################

for k = 1:samples
    # Nominal cases
    # --------------------------------------------------------------------------
    for rot_seq in rot_seq_array
        # Sample three angles form a uniform distribution [-pi,pi].
        θx = -pi + 2*pi*rand()
        θy = -pi + 2*pi*rand()
        θz = -pi + 2*pi*rand()
        Θ  = EulerAngles(θx, θy, θz, rot_seq)

        # Get the error matrix related to the conversion from DCM => Euler
        # Angles => DCM.
        error1 = angle_to_rot(DCM, dcm_to_angle(angle_to_rot(Θ),rot_seq)) -
                 angle_to_rot(DCM, θx, θy, θz, rot_seq)
        error2 = angle_to_rot(dcm_to_angle(angle_to_rot(θx, θy, θz, rot_seq),rot_seq)) -
                 angle_to_rot(Θ)

        # If everything is fine, the norm of the matrix error should be small.
        @test norm(error1) < 1e-10
        @test norm(error2) < 1e-10
        @test error1 ≈ error2
    end

    # Gimbal-lock and special cases
    # --------------------------------------------------------------------------

    # Rotation about the three axes
    # -----------------------------
    for rot_seq in [:ZYX, :ZXY, :YZX, :YXZ, :XZY, :XYZ]
        # A 90 deg rotation about the second axis yields a gimbal-lock.
        θ₁  = -π + 2π*rand()
        θ₃  = -π + 2π*rand()

        Dr1 = angle_to_dcm(θ₁, +π/2, θ₃, rot_seq)
        Dr2 = angle_to_dcm(θ₁, -π/2, θ₃, rot_seq)

        Θr1 = dcm_to_angle(Dr1, rot_seq)
        Θr2 = dcm_to_angle(Dr2, rot_seq)

        θs1, θs2 = (rot_seq in [:ZYX, :YXZ, :XZY]) ? (θ₁ - θ₃, θ₁ + θ₃) :
                                                     (θ₁ + θ₃, θ₁ - θ₃)

        # The angles must be between [-π,+π].
        θs1 = mod(θs1, 2π)
        θs2 = mod(θs2, 2π)

        (θs1 > π) && (θs1 -= 2π)
        (θs2 > π) && (θs2 -= 2π)

        @test Θr1.a1 ≈ θs1   atol = 1e-14
        @test Θr1.a2 ≈ +π/2  atol = 1e-14
        @test Θr1.a3 ≈ 0     atol = 1e-14

        @test Θr2.a1 ≈ θs2   atol = 1e-14
        @test Θr2.a2 ≈ -π/2  atol = 1e-14
        @test Θr2.a3 ≈ 0     atol = 1e-14
    end
    break

    # Rotation about two axes
    # -----------------------
    Rx = create_rotation_matrix(π, :X)
    Ry = create_rotation_matrix(π, :Y)
    Rz = create_rotation_matrix(π, :Z)
    α  = -pi + 2*pi*rand()

    # Two rotations about X
    D  = create_rotation_matrix(α, :X)
    Θ₁ = dcm_to_angle(D, :XYX)
    Θ₂ = dcm_to_angle(D, :XZX)
    @test Θ₁.a1 ≈ α
    @test Θ₁.a2 ≈ 0
    @test Θ₁.a3 ≈ 0
    @test Θ₂.a1 ≈ α
    @test Θ₂.a2 ≈ 0
    @test Θ₂.a3 ≈ 0

    Θ₁ = dcm_to_angle(D*Ry, :XYX)
    Θ₂ = dcm_to_angle(D*Rz, :XZX)
    @test Θ₁.a1 ≈ -α
    @test Θ₁.a2 ≈ π
    @test Θ₁.a3 ≈ 0
    @test Θ₂.a1 ≈ -α
    @test Θ₂.a2 ≈ π
    @test Θ₂.a3 ≈ 0

    # Two rotations about Y
    D  = create_rotation_matrix(α, :Y)
    Θ₁ = dcm_to_angle(D, :YXY)
    Θ₂ = dcm_to_angle(D, :YZY)
    @test Θ₁.a1 ≈ α
    @test Θ₁.a2 ≈ 0
    @test Θ₁.a3 ≈ 0
    @test Θ₂.a1 ≈ α
    @test Θ₂.a2 ≈ 0
    @test Θ₂.a3 ≈ 0

    Θ₁ = dcm_to_angle(D*Rx, :YXY)
    Θ₂ = dcm_to_angle(D*Rz, :YZY)
    @test Θ₁.a1 ≈ -α
    @test Θ₁.a2 ≈ π
    @test Θ₁.a3 ≈ 0
    @test Θ₂.a1 ≈ -α
    @test Θ₂.a2 ≈ π
    @test Θ₂.a3 ≈ 0

    # Two rotations about Z
    D  = create_rotation_matrix(α, :Z)
    Θ₁ = dcm_to_angle(D, :ZXZ)
    Θ₂ = dcm_to_angle(D, :ZYZ)
    @test Θ₁.a1 ≈ α
    @test Θ₁.a2 ≈ 0
    @test Θ₁.a3 ≈ 0
    @test Θ₂.a1 ≈ α
    @test Θ₂.a2 ≈ 0
    @test Θ₂.a3 ≈ 0

    Θ₁ = dcm_to_angle(D*Rx, :ZXZ)
    Θ₂ = dcm_to_angle(D*Ry, :ZYZ)
    @test Θ₁.a1 ≈ -α
    @test Θ₁.a2 ≈ π
    @test Θ₁.a3 ≈ 0
    @test Θ₂.a1 ≈ -α
    @test Θ₂.a2 ≈ π
    @test Θ₂.a3 ≈ 0
end

for k = 1:samples
    # Sample three angles form a uniform distribution [-0.0001,0.0001].
    θx = -0.0001 + 0.0002*pi*rand()
    θy = -0.0001 + 0.0002*pi*rand()
    θz = -0.0001 + 0.0002*pi*rand()
    Θ  = EulerAngles(θx, θy, θz, :XYZ)

    # Get the error between the exact rotation and the small angle
    # approximation.
    error1 = angle_to_rot(Θ) - smallangle_to_rot(DCM, Θ.a1, Θ.a2, Θ.a3)
    error2 = angle_to_rot(θx, θy, θz, :XYZ) - smallangle_to_rot(Θ.a1, Θ.a2, Θ.a3)

    # If everything is fine, the norm of the matrix error should be small.
    @test norm(error1) < 1e-9
    @test norm(error2) < 1e-9
    @test error1 ≈ error2
end

# Test orthonormalization.
let
    θx = -0.01 + 0.02*pi*rand()
    θy = -0.01 + 0.02*pi*rand()
    θz = -0.01 + 0.02*pi*rand()

    D1 = smallangle_to_rot(θx, θy, θz)
    D2 = smallangle_to_rot(θx, θy, θz; normalize = false)

    @test norm(I - D1*D1') ≈    0 atol = 1e-15
    @test norm(I - D2*D2') > 1e-7
end

# Test exceptions.
@test_throws ArgumentError angle_to_dcm(0,0,0,:xyz)
@test_throws ArgumentError angle_to_dcm(0,0,0,:zyx)
@test_throws ArgumentError angle_to_dcm(0,0,0,:xyx)
@test_throws ArgumentError angle_to_dcm(0,0,0,:abc)
@test_throws ArgumentError angle_to_quat(0,0,0,:xyz)
@test_throws ArgumentError angle_to_quat(0,0,0,:zyx)
@test_throws ArgumentError angle_to_quat(0,0,0,:xyx)
@test_throws ArgumentError angle_to_quat(0,0,0,:abc)


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

import Test
Test.@testset "Package UpROOT" begin

include("test_pycall.jl")
include("test_testfiles.jl")

end # testset


text:
using SafeTestsets

@safetestset "PolynomialChaosSurrogates" begin

using Surrogates
using PolyChaos
using Surrogates: sample, SobolSample
using SurrogatesPolyChaos
using Zygote

#1D
n = 20
lb = 0.0
ub = 4.0
f = x-> 2*x
x = sample(n,lb,ub,SobolSample())
y = f.(x)

my_pce = PolynomialChaosSurrogate(x,y,lb,ub)
val = my_pce(2.0)
add_point!(my_pce,3.0,6.0)
my_pce_changed = PolynomialChaosSurrogate(x,y,lb,ub,op=Uniform01OrthoPoly(1))

#ND
n = 60
lb = [0.0,0.0]
ub = [5.0,5.0]
f = x-> x[1]*x[2]
x = sample(n,lb,ub,SobolSample())
y = f.(x)

my_pce = PolynomialChaosSurrogate(x,y,lb,ub)
val = my_pce((2.0,2.0))
add_point!(my_pce,(2.0,3.0),6.0)

op1 = Uniform01OrthoPoly(1)
op2 = Beta01OrthoPoly(2,2,1.2)
ops = [op1,op2]
multi_poly = MultiOrthoPoly(ops,min(1,2))
my_pce_changed = PolynomialChaosSurrogate(x,y,lb,ub,op = multi_poly)

# Surrogate optimization test
lb = 0.0
ub = 15.0
p = 1.99
a = 2
b = 6
objective_function = x -> 2*x+1
x = sample(20,lb,ub,SobolSample())
y = objective_function.(x)
my_poly1d = PolynomialChaosSurrogate(x,y,lb,ub)
surrogate_optimize(objective_function,SRBF(),a,b,my_poly1d,LowDiscrepancySample(2))

lb = [0.0,0.0]
ub = [10.0,10.0]
obj_ND = x -> log(x[1])*exp(x[2])
x = sample(40,lb,ub,UniformSample())
y = obj_ND.(x)
my_polyND = PolynomialChaosSurrogate(x,y,lb,ub)
surrogate_optimize(obj_ND,SRBF(),lb,ub,my_polyND,SobolSample(),maxiters=15)

# AD Compatibility

lb = 0.0
ub = 3.0
n = 10
x = sample(n,lb,ub,SobolSample())
f = x -> x^2
y = f.(x)

# #Polynomialchaos
@testset "Polynomial Chaos" begin
	f = x -> x^2
    n = 50
    x = sample(n,lb,ub,SobolSample())
    y = f.(x)
    my_poli = PolynomialChaosSurrogate(x,y,lb,ub)
    g = x -> my_poli'(x)
    g(3.0)
end

# #PolynomialChaos
@testset "Polynomial Chaos ND" begin
    n = 50
    lb = [0.0,0.0]
    ub = [10.0,10.0]
    x = sample(n,lb,ub,SobolSample())
    f = x -> x[1]*x[2]
    y = f.(x)
    my_poli_ND = PolynomialChaosSurrogate(x,y,lb,ub)
    g = x -> Zygote.gradient(my_poli_ND,x)
    g((1.0,1.0))

    n = 10
    d = 2
    lb = [0.0,0.0]
    ub = [5.0,5.0]
    x = sample(n,lb,ub,SobolSample())
    f = x -> x[1]^2 + x[2]^2
    y1 = f.(x)
    grad1 = x -> 2*x[1]
    grad2 = x -> 2*x[2]
    function create_grads(n,d,grad1,grad2,y)
        c = 0
        y2 = zeros(eltype(y[1]),n*d)
        for i = 1:n
            y2[i+c] = grad1(x[i])
            y2[i+c+1] = grad2(x[i])
            c = c+1
        end
        return y2
    end
    y2 = create_grads(n,d,grad1,grad2,y)
    y = vcat(y1,y2)
    my_gek_ND = GEK(x,y,lb,ub)
    g = x -> Zygote.gradient(my_gek_ND,x)
    g((2.0,5.0))
end

end


text:
S = randSeisData(3, s=1.0, nx=1024)
@test mkchans(1, S) == [1]
@test mkchans(1:2, S) == [1, 2]
@test mkchans([2,3], S) == [2, 3]

S.x[2] = Float64[]
@test mkchans([2,3], S) == [3]

S.t[3] = Array{Int64, 2}(undef, 0, 2)
@test mkchans(1:3, S, f=:t) == [1,2]
@test mkchans(1:3, S) == [1,3]

S.t[3] = [1 S.t[1][1,2]; length(S.x[3]) 0]
push!(S, randSeisChannel(c=true))
@test mkchans(1:4, S, keepirr=false) == [1, 3]
@test mkchans([2,3,4], S, keepirr=false) == [3]

S = randSeisData(3, s=1.0, nx=1024)
@test get_seis_channels(S, chans=3) == [3]
@test get_seis_channels(S, chans=1:2) == [1,2]
@test get_seis_channels(S, chans=[1,3]) == [1,3]

cr = [1,2,3]
c0 = deepcopy(cr)
filt_seis_chans!(cr, S)
@test cr == c0

S.id[1] = "...0"
filt_seis_chans!(cr, S)
@test cr != c0
@test cr == [2,3]


text:
# # West Antarctica
#
# This is our first real world example. We produce a map of approximate ice velocity in West Antarctica.
# This example demonstrates:
#   * the capability of WAVI.jl in real world examples
#   * how to specify the solution mask `h_mask`, which defines which grid points are part of the ice domain.
#   * how to specify a spatially variable initial ice viscosity

# **NB** the files included here are intended as a full demonstration of WAVI.jl capability, but should not be interpreted as "ready to go" for scientific inquiry; please get in touch if you would like to use these files for science!

# ## Install dependencies
#
# First let's make sure we have all required packages installed.
# As well as WAVI and Plots for plotting, we're going to use the Downloads package to pull some data from a Github repository.

#using Pkg
#Pkg.add("https://github.com/RJArthern/WAVI.jl"), Pkg.add(Plots); Pkg.add("Downloads")
using WAVI, Plots, Downloads

# ## Reading in data
# First, let's define our grid sizes and origin. We have a 268 x 315 with, with 3km resolution.
# We're going to use 12 levels in the vertical: even though WAVI.jl is designed for the solution of depth integrated equations, it retains some information about the vertical direction e.g. in the calculation of the ice viscosity; the keyword argument `nσ`, which is passed to a `Grid` object, specifies the number of levels in the vertical. We refer to the extension of the 2D (horizontal) grid to n$\sigma$ levels in the vertical as the 3D grid, which has size nx $\times$ ny $\times$ nz.
nx = 268;
ny = 315;
nσ = 12;
x0 = -1792500.0;
y0 = -838500.0;
dx = 3000.0;
dy = 3000.0;

# We're going to read in the following data files:
#   * h_mask     : array of zeros and ones defining the ice domain (zero corresponds to out of domain, one to in domain).
#   * u_iszero   : location of grid points with zero velocity in x
#   * v_iszero   : location of grid points with zero velocity in y
#   * bed        : the bed elevation, which is a processed form of data from Bedmap 2.
#   * h          : ice thickness , defined on the 2D grid
#   * viscosity  : the initial ice viscosity, defined on the 3D grid

# We need the first three of these before we can build a grid
h_mask=Array{Float64}(undef,nx,ny);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_h_mask_clip.bin"),h_mask);
hm = ntoh.(h_mask);
hm = map.(Bool, round.(Int, hm));

u_iszero=Array{Float64}(undef,nx+1,ny);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_uiszero_clip.bin"),u_iszero);
u_iszero.=ntoh.(u_iszero);
u_iszero = map.(Bool, round.(Int, u_iszero));

v_iszero=Array{Float64}(undef,nx,ny+1);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_viszero_clip.bin"),v_iszero);
v_iszero.=ntoh.(v_iszero);
v_iszero = map.(Bool, round.(Int, v_iszero));

grid = Grid(nx = nx, ny = ny, nσ = nσ, x0 = x0, y0 = y0, dx = dx, dy = dy, h_mask = hm, u_iszero = u_iszero, v_iszero = v_iszero);

# Before moving on to the rest, let's have a look at h_mask using the `spy` function, which indicates which entries of a matrix are non-zero:
#Plots.spy(hm)
# We see a rough outline of West Antarctica! This is out ice domain. Note that indices with h_mask[i,j] = 0 are simply placeholders and do not have any physical meaning.

# Next up: the bed, which will be passed to a model via the `bed_elevation` keyword article, as usual.
bed=Array{Float64}(undef,nx,ny);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_bed_clip_noNan.bin"),bed);
bed.=ntoh.(bed)

# Let's take a look at the bed
plt = Plots.heatmap(grid.xxh[:,1]/1e3, grid.yyh[1,:]/1e3, bed',
                    xlabel = "x (km)",
                    ylabel = "y (km)",
                    colorbar_title = "bed elevation (m)",
                    title = "West Antarctica bed elevation",
                    framestyle = "box")

# Temperature, damage, viscosity, and thickness will be passed to a model via initial conditions (these quantities are time dependent)
temp=Array{Float64}(undef,nx,ny,nσ);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_3Dtemp_clip_noNan.bin"),temp);
temp.=ntoh.(temp)

damage=Array{Float64}(undef,nx,ny,nσ);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_damage3D_clip_noNan.bin"),damage);
damage.=ntoh.(damage)

h=Array{Float64}(undef,nx,ny);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_thickness_clip_noNan.bin"),h);
h.=ntoh.(h)

viscosity=Array{Float64}(undef,nx,ny,nσ);
read!(Downloads.download("https://github.com/alextbradley/WAVI_example_data/raw/main/WAIS/Inverse_3km_viscosity3D_clip_noNan.bin"),viscosity);
viscosity.=ntoh.(viscosity);

initial_conditions = InitialConditions(initial_thickness = h,initial_viscosity = viscosity,initial_temperature = temp,initial_damage = damage);

# ## Ice Velocity
# Now we're ready to make our model, which we can then use to determine the ice velocity. All physical and solver parameters take their default values
model = Model(grid = grid, bed_elevation = bed,initial_conditions= initial_conditions);

# We use the `update_state!` method to bring fields (including velocity) in line with the ice thickness:
update_state!(model);

# Now we can visualize the ice velocity:
plt = Plots.heatmap(grid.xxh[:,1]/1e3, grid.yyh[1,:]/1e3, model.fields.gh.av_speed',
                    xlabel = "x (km)",
                    ylabel = "y (km)",
                    colorbar_title = "ice speed (m/yr)",
                    title = "West Antarctica ice speed",
                    framestyle = "box")

# Ice velocities reach a maximum of approx 5km/yr on ice shelves, but are much smaller inland.


text:
"""
`find_all_combinations(vectors)`

Find all possible combinations among the vectors of values.

**Arguments**
- `vectors`: Vector of vectors representing different value sets.

**Returns**
- `combinations`: Matrix containing all combinations of values from the
    individual vectors. Each row represents one combination.

**Example**
```
julia> vectors = [[1, 2], [3, 4, 5], [6]];
julia> combinations = find_all_combinations(vectors);
julia> print(combinations);
[1 3 6; 1 4 6; 1 5 6; 2 3 6; 2 4 6; 2 5 6]
```
"""
function find_all_combinations(vectors)

    # Number of input vectors
    n = length(vectors)

    # The first column of the input matrix is the last input vector
    combinations = hcat(vectors[n])

    # Iterate the remaining vectors
    for i = (n-1):-1:1

        # Get the current vector and its length
        x = vectors[i]
        nx = length(x)

        # Get the current size of the combinations matrix
        nc = size(combinations, 1)

        # Repeat the current combinations for each element in x
        combinations = repeat(combinations, nx)

        # Add an element from x to each repetition of combinations
        x = vcat([fill(xi, nc) for xi in x]...)
        combinations = hcat(x, combinations)

    end

    return combinations

end

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
PROJECT_DIR = pwd()*"\\.."

using DataFrames
using Dates
#using HDF5, JLD

include("Params.jl")
include("Solution.jl")
include("LocalSearch.jl")
include("GeneticAlgorithm.jl")

#allinstances = readdir("../../instances/");
#allsolutions = readdir("../../solutions/");

#println("Number of instance files: ",length(allinstances))
#println("Number of solution files: ",length(allsolutions))

df = DataFrame(
    instance = String[],
    size = Int64[],
    seed = Int64[],
    finalCost = Int64[],
    populationSize = Int64[],
    maxGenerations = Int64[],
    recombineOp = Function[],
    distanceFunc = Function[],
    entropyReduceFunc = Function[],
    mutate = Any[],
    selection = Any[],
    executionTime = Float64[],
    numIndividuals = Int64[],
    finalEntropy = Float64[],
    numGenerations = Int64[],
    bestKnownCost = Any[],
    gap = Any[]
)

timenow = Dates.now()
start_datetime = Dates.format(timenow, "yyyy-mm-dd_HHMMSS")
experiment_dir = joinpath(PROJECT_DIR,"experiments",start_datetime)
mkdir(experiment_dir)

results_file = "aggregate_results.txt"

cpu_time = 120

cols_tuple = Tuple(string(sym) for sym in describe(df).variable)
open(joinpath(experiment_dir, results_file),"a+") do txtfile
   write(txtfile,join(cols_tuple,"\t"),"\n")
end

# Grid search of parameters
n_seed = 1:5
n_maxGenerations = 10^10
n_populationSize = [50]
n_recombineOp = [order1cx,PMXcx,cyclecx]
n_distanceFunc = [hamming]
n_entropyReduceFunc = [mean]
n_mutate = [cyclicHeuristic]
n_selection = [fitnessSelection,biasedFitnessSelection]

allinstances = [
    "wil100.dat","tho40.dat","ste36c.dat","tai60b.dat","tai80b.dat",
    "sko49.dat","sko100d.dat","nug30.dat","tai150b.dat",#"tho150.dat",
    #"chr20c.dat","wil50.dat","lipa40b.dat","lipa50a.dat","lipa90b.dat"
    #"els19.dat"
    #"bur26a.dat","bur26e.dat","esc16h.dat","esc32e.dat","esc128.dat","chr15a.dat","nug24.dat","scr20.dat"
]

for gaParams in Iterators.product(
        n_seed,n_maxGenerations,n_populationSize,n_recombineOp,n_distanceFunc,
        n_entropyReduceFunc,n_mutate,n_selection,allinstances)
	# Unpack parameters for GA
    seed,maxGenerations,populationSize,recombine,distFunc,
        entropyReduce, mutate, selection, instance_name = gaParams

    println(gaParams)
    println("$instance_name; seed=$seed; maxGen=$maxGenerations; popSize=$populationSize")
    instance_path = joinpath(PROJECT_DIR, "instances",instance_name)
    output_path = joinpath(experiment_dir,split(instance_name,'.')[1]*".jld")
    params = Params(instance_path, output_path, seed, cpu_time)

    Random.seed!(seed)
    solGA, results = runGA(params, maxGenerations, populationSize,
                           recombineOp=recombine,
                           distanceFunc=distFunc,
                           entropyReduceFunc=entropyReduce,
                           mutate=mutate,
                           selectionOp=selection,
                           selectionOpts=Dict{String,Any}("fitnessThreshold"=>0.45))

    println("GA solution:")
    printSolution(solGA)
    println("Gap=",results["gap"][end])
    println()

    if occursin(" ", params.pathToOutput)
        throw("Your project path must not contain any spaces, because of the behaviour of the JLD package.")
    end
    #save(params.pathToOutput, "data", results) # for some unknown reason this does not work
    #save(joinpath("..\\..","experiments",start_datetime,split(instance_name,'.')[1]*".jld"), "data", results)

    new_tuple = (instance_name,
        params.datasetSize,
        seed,
        results["cost"][end],
        results["populationSize"],
        results["maxGenerations"],
        recombine,
        distFunc,
        entropyReduce,
        mutate,
        selection,
        results["time"],
        results["numIndividuals"][end],
        results["entropy"][end],
        results["generations"],
        params.solutionCost,
        results["gap"][end])
    push!(df, new_tuple)

    open(joinpath(experiment_dir, results_file),"a+") do txtfile
       write(txtfile,join(new_tuple,"\t"),"\n")
    end
end

println("Finished execution.")

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics


text:
module TestReferenceables

using Dagger
using FLoops
using FoldsDagger
using Referenceables
using Test

function test_in_place_module()
    xs = distribute(collect(1:4), Blocks(2))
    @floop DaggerEx() for x in referenceable(xs)
        x[] += 1
    end
    @test collect(xs) == 2:5
end

function test_copy()
    xs = distribute(1:4, Blocks(2))
    ys = distribute(zeros(4), Blocks(2))
    @floop DaggerEx() for (y, x) in zip(referenceable(ys), xs)
        y[] = x
    end
    @test collect(ys) == 1:4
end

function test_fan_out()
    xs = distribute(1:4, Blocks(2))
    ys = distribute(zeros(4), Blocks(2))
    zs = distribute(zeros(4), Blocks(2))
    @floop DaggerEx() for (z, y, x) in zip(referenceable(zs), referenceable(ys), xs)
        y[] = x
        z[] = 2x
    end
    @test collect(ys) == 1:4
    @test collect(zs) == 2:2:8
end

end


text:
"""
    unsafe_createstylemanager(styletable = GDALStyleTable(C_NULL))

OGRStyleMgr factory.

### Parameters
* `styletable`: OGRStyleTable or NULL if not working with a style table.

### Returns
an handle to the new style manager object.
"""
unsafe_createstylemanager(styletable = GDALStyleTable(C_NULL)) =
    StyleManager(GDAL.ogr_sm_create(styletable))

"""
Destroy Style Manager.

### Parameters
* `stylemanager`: handle to the style manager to destroy.
"""
function destroy(sm::StyleManager)
    GDAL.ogr_sm_destroy(sm.ptr)
    sm.ptr = C_NULL
end

"""
    initialize!(stylemanager::StyleManager, feature::Feature)

Initialize style manager from the style string of a feature.

### Parameters
* `stylemanager`: handle to the style manager.
* `feature`: handle to the new feature from which to read the style.

### Returns
the style string read from the feature, or NULL in case of error.
"""
initialize!(stylemanager::StyleManager, feature::Feature) =
    GDAL.ogr_sm_initfromfeature(stylemanager.ptr, feature.ptr)

"""
    initialize!(stylemanager::StyleManager, stylestring = C_NULL)

Initialize style manager from the style string.

### Parameters
* `stylemanager`: handle to the style manager.
* `stylestring`: the style string to use (can be NULL).
### Returns
`true` on success, `false` on errors.
"""
initialize!(stylemanager::StyleManager, stylestring = C_NULL) =
    Bool(GDAL.ogr_sm_initstylestring(stylemanager.ptr, stylestring))

"""
    npart(stylemanager::StyleManager)
    npart(stylemanager::StyleManager, stylestring::AbstractString)

Get the number of parts in a style.

### Parameters
* `stylemanager`: handle to the style manager.
* `stylestring`: (optional) the style string on which to operate. If NULL then
                 the current style string stored in the style manager is used.

### Returns
the number of parts (style tools) in the style.
"""
function npart end

npart(stylemanager::StyleManager) =
    GDAL.ogr_sm_getpartcount(stylemanager.ptr, C_NULL)

npart(stylemanager::StyleManager, stylestring::AbstractString) =
    GDAL.ogr_sm_getpartcount(stylemanager.ptr, stylestring)

"""
    unsafe_getpart(stylemanager::StyleManager, id::Integer, stylestring = C_NULL)

Fetch a part (style tool) from the current style.

### Parameters
* `stylemanager`: handle to the style manager.
* `id`: the part number (0-based index).
* `stylestring`: (optional) the style string on which to operate. If not
    provided, then the current style string stored in the style manager is used.

### Returns
OGRStyleToolH of the requested part (style tools) or NULL on error.
"""
unsafe_getpart(stylemanager::StyleManager, id::Integer, stylestring = C_NULL) =
    StyleTool(GDAL.ogr_sm_getpart(stylemanager.ptr, id, stylestring))

"""
    addpart!(stylemanager::StyleManager, styletool::StyleTool)

Add a part (style tool) to the current style.

### Parameters
* `stylemanager`: handle to the style manager.
* `styletool`: the style tool defining the part to add.

### Returns
`true` on success, `false` on errors.
"""
addpart!(stylemanager::StyleManager, styletool::StyleTool) =
    Bool(GDAL.ogr_sm_addpart(stylemanager.ptr, styletool.ptr))

"""
    addstyle!(stylemanager::StyleManager, stylename, stylestring)

Add a style to the current style table.

### Parameters
* `stylemanager`: handle to the style manager.
* `stylename`: the name of the style to add.
* `stylestring`: (optional) the style string to use, or (if not provided) to use
    the style stored in the manager.

### Returns
`true` on success, `false` on errors.
"""
function addstyle!(
        stylemanager::StyleManager,
        stylename::AbstractString,
        stylestring::AbstractString
    )
    return Bool(GDAL.ogr_sm_addstyle(stylemanager.ptr, stylename, stylestring))
end

addstyle!(stylemanager::StyleManager, stylename::AbstractString) =
    Bool(GDAL.ogr_sm_addstyle(stylemanager.ptr, stylename, C_NULL))

"""
    unsafe_createstyletool(classid::OGRSTClassId)

OGRStyleTool factory.

### Parameters
* `classid`: subclass of style tool to create. One of OGRSTCPen (1),
             OGRSTCBrush (2), OGRSTCSymbol (3) or OGRSTCLabel (4).

### Returns
an handle to the new style tool object or NULL if the creation failed.
"""
unsafe_createstyletool(classid::OGRSTClassId) =
    StyleTool(GDAL.ogr_st_create(classid))

"""
Destroy Style Tool.

### Parameters
* `styletool`: handle to the style tool to destroy.
"""
function destroy(styletool::StyleTool)
    GDAL.ogr_st_destroy(styletool.ptr)
    styletool.ptr = C_NULL
end

"""
    gettype(styletool::StyleTool)

Determine type of Style Tool.

### Parameters
* `styletool`: handle to the style tool.

### Returns
the style tool type, one of OGRSTCPen (1), OGRSTCBrush (2), OGRSTCSymbol (3) or
OGRSTCLabel (4). Returns OGRSTCNone (0) if the OGRStyleToolH is invalid.
"""
gettype(styletool::StyleTool) = GDAL.ogr_st_gettype(styletool.ptr)

"""
    getunit(styletool::StyleTool)

Get Style Tool units.

### Parameters
* `styletool`: handle to the style tool.

### Returns
the style tool units.
"""
getunit(styletool::StyleTool) = GDAL.ogr_st_getunit(styletool.ptr)

"""
    setunit!(styletool::StyleTool, newunit::OGRSTUnitId, scale::Real)

Set Style Tool units.

### Parameters
* `styletool`: handle to the style tool.
* `newunit`: the new unit.
* `scale`: ground to paper scale factor.
"""
setunit!(styletool::StyleTool, newunit::OGRSTUnitId, scale::Real) =
    GDAL.ogr_st_setunit(styletool.ptr, newunit, scale)

"""
    asstring(styletool::StyleTool, id::Integer)
    asstring(styletool::StyleTool, id::Integer, nullflag::Ref{Cint})

Get Style Tool parameter value as a string.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `nullflag`: pointer to an integer that will be set to `true` or `false` to
        indicate whether the parameter value is NULL.

### Returns
the parameter value as a string and sets `nullflag`.
"""
asstring(styletool::StyleTool, id::Integer, nullflag::Ref{Cint}) =
    GDAL.ogr_st_getparamstr(styletool.ptr, id, nullflag)

asstring(styletool::StyleTool, id::Integer) =
    asstring(styletool, id, Ref{Cint}(0))

"""
    asint(styletool::StyleTool, id::Integer, nullflag = Ref{Cint}(0))

Get Style Tool parameter value as an integer.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `nullflag`: pointer to an integer that will be set to `true` or `false` to
        indicate whether the parameter value is NULL.

### Returns
the parameter value as an integer and sets `nullflag`.
"""
asint(styletool::StyleTool, id::Integer, nullflag::Ref{Cint} = Ref{Cint}(0)) =
    GDAL.ogr_st_getparamnum(styletool.ptr, id, nullflag)

"""
    asdouble(styletool::StyleTool, id::Integer, nullflag = Ref{Cint}(0))

Get Style Tool parameter value as a double.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `nullflag`: pointer to an integer that will be set to `true` or `false` to
        indicate whether the parameter value is NULL.

### Returns
the parameter value as a double and sets `nullflag`.
"""
function asdouble(
        styletool::StyleTool,
        id::Integer,
        nullflag::Ref{Cint} = Ref{Cint}(0)
    )
    return GDAL.ogr_st_getparamdbl(styletool.ptr, id, nullflag)
end

"""
    setparam!(styletool::StyleTool, id::Integer, value)

Set Style Tool parameter value.

### Parameters
* `styletool`: handle to the style tool.
* `id`: the parameter id from the enumeration corresponding to the type of this
        style tool (one of the OGRSTPenParam, OGRSTBrushParam, OGRSTSymbolParam
        or OGRSTLabelParam enumerations)
* `value`: the new parameter value, can be an Integer, Float64, or AbstactString
"""
function setparam! end

setparam!(styletool::StyleTool, id::Integer, value::AbstractString) =
    GDAL.ogr_st_setparamstr(styletool.ptr, id, value)

setparam!(styletool::StyleTool, id::Integer, value::Integer) =
    GDAL.ogr_st_setparamnum(styletool.ptr, id, value)

setparam!(styletool::StyleTool, id::Integer, value::Float64) =
    GDAL.ogr_st_setparamdbl(styletool.ptr, id, value)

"""
    getstylestring(styletool::StyleTool)

Get the style string for this Style Tool.

### Parameters
* `styletool`: handle to the style tool.

### Returns
the style string for this style tool or "" if the styletool is invalid.
"""
getstylestring(styletool::StyleTool) = GDAL.ogr_st_getstylestring(styletool.ptr)

"""
    toRGBA(styletool::StyleTool, color::AbstractString)

Return the r,g,b,a components of a color encoded in #RRGGBB[AA] format.

### Parameters
* `styletool`: handle to the style tool.
* `pszColor`: the color to parse

### Returns
(R,G,B,A) tuple of Cints.
"""
function toRGBA(styletool::StyleTool, color::AbstractString)
    red = Ref{Cint}(0)
    green = Ref{Cint}(0)
    blue = Ref{Cint}(0)
    alpha = Ref{Cint}(0)
    result = Bool(GDAL.ogr_st_getrgbfromstring(styletool.ptr, color, red, green,
        blue, alpha))
    result || error("Error in getting RGBA from Styletool")
    return (red[], green[], blue[], alpha[])
end

"""
    unsafe_createstyletable()

OGRStyleTable factory.

### Returns
an handle to the new style table object.
"""
unsafe_createstyletable() = StyleTable(GDAL.ogr_stbl_create())

"""
Destroy Style Table.

### Parameters
* `styletable`: handle to the style table to destroy.
"""
function destroy(st::StyleTable)
    GDAL.ogr_stbl_destroy(st.ptr)
    st.ptr = C_NULL
end

"""
    addstyle!(styletable::StyleTable, stylename, stylestring)

Add a new style in the table.

### Parameters
* `styletable`: handle to the style table.
* `name`: the name the style to add.
* `stylestring`: the style string to add.

### Returns
`true` on success, `false` on error
"""
function addstyle!(
        styletable::StyleTable,
        stylename::AbstractString,
        stylestring::AbstractString
    )
    return Bool(GDAL.ogr_stbl_addstyle(styletable.ptr, stylename, stylestring))
end

"""
    savestyletable(styletable::StyleTable, filename::AbstractString)

Save a style table to a file.

### Parameters
* `styletable`: handle to the style table.
* `filename`: the name of the file to save to.

### Returns
`true` on success, `false` on error
"""
savestyletable(styletable::StyleTable, filename::AbstractString) =
    Bool(GDAL.ogr_stbl_savestyletable(styletable.ptr, filename))

"""
    loadstyletable!(styletable::StyleTable, filename::AbstractString)

Load a style table from a file.

### Parameters
* `styletable`: handle to the style table.
* `filename`: the name of the file to load from.

### Returns
`true` on success, `false` on error
"""
loadstyletable!(styletable::StyleTable, filename::AbstractString) =
    Bool(GDAL.ogr_stbl_loadstyletable(styletable.ptr, filename))

"""
    findstylestring(styletable::StyleTable, name::AbstractString)

Get a style string by name.

### Parameters
* `styletable`: handle to the style table.
* `name`: the name of the style string to find.

### Returns
the style string matching the name or NULL if not found or error.
"""
findstylestring(styletable::StyleTable, name::AbstractString) =
    GDAL.ogr_stbl_find(styletable.ptr, name)

"""
    resetreading!(styletable::StyleTable)

Reset the next style pointer to 0.

### Parameters
* `styletable`: handle to the style table.
"""
resetreading!(styletable::StyleTable) =
    GDAL.ogr_stbl_resetstylestringreading(styletable.ptr)

"""
    nextstyle(styletable::StyleTable)

Get the next style string from the table.

### Parameters
* `styletable`: handle to the style table.

### Returns
the next style string or NULL on error.
"""
nextstyle(styletable::StyleTable) = GDAL.ogr_stbl_getnextstyle(styletable.ptr)

"""
    laststyle(styletable::StyleTable)

Get the style name of the last style string fetched with OGR_STBL_GetNextStyle.

### Parameters
* `styletable`: handle to the style table.

### Returns
the Name of the last style string or NULL on error.
"""
laststyle(styletable::StyleTable) =
    GDAL.ogr_stbl_getlaststylename(styletable.ptr)


text:
export Overwrite

const OVERWRITE_TYPES = ["role", "member"]

"""
An [`Overwrite`](@ref)'s type.
More details [here](https://discordapp.com/developers/docs/resources/audit-log#audit-log-entry-object-optional-audit-entry-info).
"""
@enum OverwriteType OT_ROLE OT_MEMBER OT_UNKNOWN
@boilerplate OverwriteType :export

function OverwriteType(x::AbstractString)
    i = findfirst(isequal(x), OVERWRITE_TYPES)
    return i === nothing ? OT_UNKNOWN : OverwriteType(i - 1)
end

Base.string(x::OverwriteType) = x === OT_UNKNOWN ? "" : OVERWRITE_TYPES[Int(x) + 1]
JSON.lower(x::OverwriteType) = string(x)

"""
A permission overwrite.
More details [here](https://discordapp.com/developers/docs/resources/channel#overwrite-object).
"""
struct Overwrite
    id::Snowflake
    type::OverwriteType
    allow::Int
    deny::Int
end
@boilerplate Overwrite :constructors :docs :lower :merge :mock


text:
module HypergraphModularity

import Combinatorics
import Distributions
import LinearAlgebra
import Parameters
import Random
import SparseArrays
import StatsBase
import Base
import NLopt
import Optim
import SpecialFunctions
import DelimitedFiles  #

# eventually move this and the contents of the read_data
# module over to a separate script, probably in test,
# as it is not part of the core hypergraph modularity
# suite

include("omega.jl")
include("HSBM.jl")

include("vol.jl")

include("utils.jl")

include("dict_ops.jl")

include("inference.jl")
include("objectives.jl")

include("graph_louvain.jl")
include("louvain_utils.jl")
include("general_hypergraph_louvain.jl")
include("read_data.jl")          # we should move this out of the package
include("analysis_helpers.jl")   # this too

include("test_funs.jl")

include("AON_hyperlouvain.jl")

export partitionize
export sampleSBM

export counting_coefficient

export CliqueExpansion
export CliqueExpansionModularity
export HyperLouvain
export SuperNodeLouvain
export SuperNodeStep

export hypergraph

export buildΩ

export evalSums
export evalCuts
export evalConstants
export increments
export addIncremens

export logLikelihood
export modularity
export formObjective

export coordinateAscent

export estimateΩEmpirically

export read_hypergraph_data
export read_hypergraph_label_names
export read_hypergraph_labels
export read_hypergraph_edges

export L

export downSampleEdges!
export subhypergraph

export computeDyadicResolutionParameter
export dyadicModularity
export dyadicLogLikelihood

export partitionIntensityFunction
export IntensityFunction
export allOrNothingIntensityFunction
export sumOfExteriorDegreesIntensityFunction

export learnParameters

export AON_Inputs
export SuperNode_PPLouvain

export mutualInformation

export hyperedges
export read_stats_data
export read_cora
export readTemporalData
export subHypergraph
export dyadicLogLikelihoodDirect

export projectedGraph
export kcore
export removeEdges!
end # module

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: PhilChodrow/HypergraphModularity

max_forks_repo_name: PhilChodrow/HypergraphModularity


text:
#---------------------------------------------------------------------#
#This routine advances the solution in time using
#a simple General Order RK method.
#Written by F.X. Giraldo / P.R. Mugg on 9/21/19
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
function WS_general_rk(q0,c,I,Δt,time)

    #build RK Coefficients
    αi = zeros(I+1,1)
    for i = 1:I
        αi[i+1] = 1.0/(I-i+1)
    end

    qi = q0
    (rhs,f,g) = rhs_function(qi,c,time)
    for i = 1:I
        #Update
        qi = q0 + Δt*αi[i+1]*rhs

        #Construct RHS
        (rhs,f,g) = rhs_function(qi,c,(time + Δt*αi[i+1]))
    end #for i

    #final update
    return qp = qi

end #end function


text:
function hausdorff(ŷ, y, ŷ_dtm, y_dtm)
    M = (ŷ .- y) .^ 2 .* (ŷ_dtm .^ 2 .+ y_dtm .^ 2)
    return loss = mean(M)
end

# TODO: make functional tullio.jl hausdorff loss
# function hausdorff(ŷ, y, ŷ_dtm, y_dtm)
#     @tullio tot :=
#         (ŷ[i, j, k, c, b] .- y[i, j, k, c, b])^2 *
#         (ŷ_dtm[i, j, k, c, b]^2 + y_dtm[i, j, k, c, b]^2)
#     return loss = tot / length(y)
# end


text: f(x,y) = 2x+3y


text:
import MathProgBase
@testset "Linear tests" begin
    include(joinpath(dirname(dirname(pathof(MathProgBase))), "test", "linproginterface.jl"))
    linprogsolvertest(SDPASolver(), 1e-5)
end
@testset "Conic tests" begin
    include(joinpath(dirname(dirname(pathof(MathProgBase))), "test", "conicinterface.jl"))

    @testset "Conic linear tests" begin
        coniclineartest(SDPASolver(), duals=true, tol=1e-5)
    end

    @testset "Conic SOC tests" begin
        conicSOCtest(SDPASolver(), duals=true, tol=1e-5)
    end

    @testset "Conic SOC rotated tests" begin
        conicSOCRotatedtest(SDPASolver(), duals=true, tol=1e-5)
    end

    @testset "Conic SDP tests" begin
        conicSDPtest(SDPASolver(), duals=false, tol=1e-5)
    end
end

using SemidefiniteModels
@testset "MPB interface" begin
    solver = SDPASolver()
    @test supportedcones(solver) == [:Free,:Zero,:NonNeg,:NonPos,:SOC,:RSOC,:SDP]
    m = SDModel(solver)
    @test_throws ErrorException loadproblem!(m, "in.dat-s")
    # If I put 0 as third argument, the .cov coverage files get deleted
    # Weirdest thing ever :o
    loadproblem!(m, [1], 1)
    @test_throws ErrorException setvartype!(m, :Int, 1, 1, 1)
    @test status(m) == :Uninitialized
end


text:
module Scenes

using VsRecorderBase
using ..VsRecorder
using ..VsRecorder: AbstractPokemonScene,
    PokemonContext, BattleContext,
    update_team!,
    default_context,
    ParseType,
    PARSE_NONE, PARSE_MINIMAL, PARSE_OPPONENT, PARSE_BOTH_PLAYERS
using ..Data

include("common.jl")
include("searching.jl")
include("team_preview.jl")
include("battle_main.jl")
include("team_view.jl")

const AvailableScenes = Type[
    Searching,
    TeamPreview, TeamPreviewSelecting,
    BattleMain,
    TeamView,
]

end


text:
using Functors, Test
using Zygote

@testset "Functors.jl" begin

  include("basics.jl")
  include("update.jl")
end


text:
using Test, Dates, TickTock

tick()

sleep(1)

@test typeof(peektimer()) == Float64

@async alarm(now() + Dates.Minute(0) + Dates.Second(10),
    action=() ->
        begin
            println("test alarm has fired")
        end)

alarmmessage() = println("This is an alarm. Do not be alarmed.")

@async alarm(now() + Dates.Second(5), action=alarmmessage)

alarmlist()

@test sprint(show, laptimer(), context=:compact => true) == "nothing"

@test tok() > 1.0

tick()

@test sprint(show, tock(), context=:compact => true) == "nothing"

println("Make 10 timers")

for i in 1:10
    tick()
    sleep(1)
end

@test length(TickTock.gettimers()) == 10

println("Check that there are 10 timers")

println("wait a second")

sleep(1)

println("Print all the timers")

laptimer()

for i in 1:10
    peektimer()
end

@show peektimers()

println("Check that the most recent timer is more than 1 second")
println("  ", peektimer())
@test peektimer() > 1.0

println("Finish 10 timers and show canonical")

for i in 1:10
    tock()
end

@test_throws ErrorException tok()

max_stars_repo_name: cormullion/TickTock

max_issues_repo_name: cormullion/TickTock

max_forks_repo_name: cormullion/TickTock


text:
using ProgressMeter

T = 50
N = 50

pomdp = BabyPOMDP()

# bds = IndependentBounds(PORollout(FeedWhenCrying(), PreviousObservationUpdater()), 0.0)
bds = IndependentBounds(SemiPORollout(FeedWhenCrying()), 0.0)

solver = AdaOPSSolver(epsilon_0=0.1,
                      m_min=200,
                      max_depth=50,
                      bounds=bds,
                      T_max=Inf,
                      max_trials=500,
                      rng=MersenneTwister(4)
                     )

rsum = 0.0
fwc_rsum = 0.0
@showprogress for i in 1:N
    planner = solve(solver, pomdp)
    sim = RolloutSimulator(max_steps=T, rng=MersenneTwister(i))
    fwc_sim = deepcopy(sim)
    global rsum += simulate(sim, pomdp, planner)
    global fwc_rsum += simulate(fwc_sim, pomdp, FeedWhenCrying())
end

@show rsum/N
@show fwc_rsum/N


text:
using JLD2, Printf

@load "MNIST_1000.jld2" train_x train_y

file = open("data.txt", "a")

print(file, "data = [")

for i in 0:9
    count = 0
    for j in 1:1000
        if train_y[j] == i
            data = reshape(train_x[:,:,j], 784)
            print(file, "[")
            for k in 1:784
                @printf(file, "%f, ", data[k])
            end
            print(file, "],\n")
            count += 1
        end
        if count == 10
            break
        end
    end
end

print(file, "]")

close(file)


text:
# Polynomial mean function

@doc """
# Description
Constructor for the Polynomial mean function

m(x) = Σ_d xᵈβᵀ_d
# Arguments:
* `β::Matrix{Float64}`: A Dxd matrix of coefficients where D is the dimension of xᵈ and d is the degree of the polynomial
""" ->
type MeanPoly <: Mean
    β::Matrix{Float64}    # Polynomial coefficients
    dim::Int              # Dimension
    deg::Int              # Polynomial degree
    MeanPoly(β::Matrix{Float64}) = new(β,size(β, 1), size(β, 2))
end

function mean(mPoly::MeanPoly,x::Matrix{Float64})
    dim, nobsv = size(x)
    dim == mPoly.dim || throw(ArgumentError("Observations and mean function have inconsistent dimensions"))
    z = zeros(nobsv)
    for i in 1:nobsv
        for j in 1:mPoly.deg
            z[i] += dot(x[:, i].^j, mPoly.β[:,j])
        end
    end
    return z
end

get_params(mPoly::MeanPoly) = vec(mPoly.β)
get_param_names(mPoly::MeanPoly) = get_param_names(mPoly.β, :β)
num_params(mPoly::MeanPoly) = mPoly.dim * mPoly.deg
function set_params!(mPoly::MeanPoly, hyp::Vector{Float64})
    num_param = mPoly.dim * mPoly.deg
    length(hyp) == num_param || throw(ArgumentError("Polynomial mean function has $(num_param) parameters"))
    mPoly.β = reshape(hyp,mPoly.dim,mPoly.deg)
end

function grad_mean(mPoly::MeanPoly, x::Vector{Float64})
    dM_theta = Array(Float64,mPoly.dim,mPoly.deg)

    for i in 1:mPoly.dim
        for j in 1:mPoly.deg
            dM_theta[i,j] = x[i].^j
        end
    end
    return vec(dM_theta)
end


text:
using Nextion
using Nextion: NexID
using Test

@testset "page" begin

    nexSerial = NexSerial("/dev/ttyUSB0")

    # Nextion pages are 0-indexed
    # nexPages is a Julia Vector and is 1-indexed

    nexPages = NexPage[]
    pages_nb = 22
    for i in 0:pages_nb-1
        nexPage = NexPage(nexSerial, Name("page$i"), pid=PageID(i), cid=ComponentID(0))
        push!(nexPages, nexPage)
    end

    for i in 1:pages_nb
        show(nexPages[i], PageID)
        sleep(0.5)
    end

    sleep(2)

    nexPage = nexPages[1]  # page0
    show(nexPage)
    @test PageID(NexID(nexPage)) == PageID(0x00)

    sleep(2)

    @test PageID(nexSerial) == PageID(0x00)  # return current page (send command sendme and wait response)

    sleep(0.1)

    @test isshown(nexPages[1])  # Page 1 (ie page0) is visible
    @test !isshown(nexPages[2])  # Page 2 (ie page1) is not visible

    close(nexSerial)
    @test true

end


text:
struct Sphere{T} <: AbstractVolumePrimitive{T, 3}
    org::CartesianPoint{T}
    r::T
end

function in(pt::CartesianPoint{T}, s::Sphere{T})::Bool where {T}
    v::CartesianVector{T} = s.org - pt
    return norm(v) <= s.r
end
function in(pt::CylindricalPoint{T}, s::Sphere{T})::Bool where {T}
    return in(CartesianPoint(pt), s)
end

function Sphere{T}(dict::Union{Dict{Any, Any}, Dict{String, Any}}, inputunit_dict::Dict{String,Unitful.Units})::Sphere{T} where {T <: SSDFloat}
    translate_vector = _get_translate_vector(T, dict, inputunit_dict)
    org::CartesianPoint{T} = ismissing(translate_vector) ? CartesianPoint{T}(0, 0, 0) : translate_vector
    r::T = ustrip(uconvert(u"m", dict["r"] * inputunit_dict["length"]))
    return Sphere{T}(org, r)
end

function Geometry(T::DataType, t::Val{:sphere}, dict::Dict{Any, Any}, inputunit_dict::Dict{String,Unitful.Units})
    return Sphere{T}(dict, inputunit_dict)
end

function (+)(s::Sphere{T}, translate::Union{CartesianVector{T}, Missing})::Sphere{T} where {T <: SSDFloat}
    return ismissing(translate) ? s : Sphere{T}( s.org +  translate, s.r )
end

@recipe function f(s::Sphere{T};) where {T <: SSDFloat}
    @series begin
        pts = []
        label --> ""
        for φ in range(T(0), length = 36, stop = T(2π))
            pt = s.org + CartesianVector{T}( s.r * cos(φ), s.r * sin(φ), 0  )
            push!(pts, pt)
        end
        lines = LineSegment{T, 3, :cartesian}[]
        for i in 1:length(pts)-1
            push!(lines, LineSegment(pts[i+1], pts[i]))
        end
        lines
    end
    @series begin
        label --> ""
        pts = []
        for φ in range(T(0), length = 36, stop = T(2π))
            pt = s.org + CartesianVector{T}( s.r * cos(φ), 0, s.r * sin(φ) )
            push!(pts, pt)
        end
        lines = LineSegment{T, 3, :cartesian}[]
        for i in 1:length(pts)-1
            push!(lines, LineSegment(pts[i+1], pts[i]))
        end
        lines
    end
    @series begin
        pts = []
        label --> ""
        for φ in range(T(0), length = 36, stop = T(2π))
            pt = s.org + CartesianVector{T}( 0, s.r * cos(φ), s.r * sin(φ) )
            push!(pts, pt)
        end
        lines = LineSegment{T, 3, :cartesian}[]
        for i in 1:length(pts)-1
            push!(lines, LineSegment(pts[i+1], pts[i]))
        end
        lines
    end
end

# For proper grid creation we also need the function get_important_points:
function get_important_points(s::Sphere{T}, ::Val{:r})::Vector{T} where {T <: SSDFloat}
    v::Vector{T} = geom_round.(T[ s.org.x - s.r, s.org.x, s.org.x + s.r,
                                  s.org.y - s.r, s.org.y, s.org.y + s.r ])
    return findall(r -> r >= 0 , v)
end
function get_important_points(s::Sphere{T}, ::Val{:φ})::Vector{T} where {T <: SSDFloat}
    return T[ ]
end
function get_important_points(s::Sphere{T}, ::Val{:z})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[ s.org.z - s.r, s.org.z, s.org.z + s.r])
end
function get_important_points(s::Sphere{T}, ::Val{:x})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[ s.org.x - s.r, s.org.x, s.org.x + s.r])
end
function get_important_points(s::Sphere{T}, ::Val{:y})::Vector{T} where {T <: SSDFloat}
    return geom_round.(T[ s.org.y - s.r, s.org.y, s.org.y + s.r])
end

# and a sample function to paint the primitive on the grid (necessary if the object is small)
function sample(s::Sphere{T}, stepsize::Vector{T})::Vector{CartesianPoint{T}}  where {T <: SSDFloat}
    samples::Vector{CartesianPoint{T}} = CartesianPoint{T}[]
    xarr::Vector{T} = get_important_points(s, Val{:x}())
    yarr::Vector{T} = get_important_points(s, Val{:y}())
    zarr::Vector{T} = get_important_points(s, Val{:z}())
    for x in xarr
        for y in yarr
            for z in zarr
                pt::CartesianPoint{T} = CartesianPoint{T}(x, y, z)
                push!(samples, pt)
            end
        end
    end
    return samples
end


text:
using Test, Pkg
using MPI

if haskey(Pkg.installed(), "CuArrays")
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end
using DoubleFloats

MPI.Init()

comm = MPI.COMM_WORLD
sz = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

root = sz-1
isroot = rank == root

val = isroot ? sum(0:sz-1) : nothing
@test MPI.Reduce(rank, MPI.SUM, root, comm) == val

val = isroot ? sz-1 : nothing
@test MPI.Reduce(rank, MPI.MAX, root, comm) == val

val = isroot ? 0 : nothing
@test MPI.Reduce(rank, MPI.MIN, root, comm) == val

val = isroot ? sz : nothing
@test MPI.Reduce(1, +, root, comm) == val

mesg = ArrayType(1.0:5.0)
sum_mesg = MPI.Reduce(mesg, +, root, comm)
if isroot
    @test sum_mesg isa ArrayType{Float64}
    @test sum_mesg == sz .* mesg
end

if ArrayType != Array || Sys.iswindows() && Sys.WORD_SIZE == 32
    operators = [MPI.SUM, +]
else
    operators = [MPI.SUM, +, (x,y) -> 2x+y-x]
end

for T = [Int]
    for dims = [1, 2, 3]
        send_arr = ArrayType(zeros(T, Tuple(3 for i in 1:dims)))
        send_arr[:] .= 1:length(send_arr)

        for op in operators

            # Non allocating version
            recv_arr = ArrayType{T}(undef, size(send_arr))
            MPI.Reduce!(send_arr, recv_arr, op, root, MPI.COMM_WORLD)
            if isroot
                @test recv_arr == sz .* send_arr
            end

            # Assertions when output buffer too small
            recv_arr = ArrayType{T}(undef, size(send_arr).-1)
            if isroot
                @test_throws AssertionError MPI.Reduce!(send_arr, recv_arr,
                                                        length(send_arr),
                                                        op, root,
                                                        MPI.COMM_WORLD)
            end

            # IN_PLACE
            recv_arr = copy(send_arr)
            MPI.Reduce_in_place!(recv_arr, length(recv_arr), op, root, MPI.COMM_WORLD)
            if isroot
                @test recv_arr == sz .* send_arr
            end

            # Allocating version
            val = MPI.Reduce(2, op, root, MPI.COMM_WORLD)
            if isroot
                @test val == sz*2
            end

            recv_arr = MPI.Reduce(send_arr, op, root, MPI.COMM_WORLD)
            if isroot
                @test recv_arr isa ArrayType{T}
                @test recv_arr == sz .* send_arr
            end

            # Allocating, Subarray
            recv_arr = MPI.Reduce(view(send_arr, 2:3), op, root, MPI.COMM_WORLD)
            if isroot
                @test recv_arr isa ArrayType{T}
                @test recv_arr == sz .* view(send_arr, 2:3)
            end
        end
    end
end

MPI.Barrier( MPI.COMM_WORLD )

if !(Sys.iswindows() && Sys.WORD_SIZE == 32)

    send_arr = [Double64(i)/10 for i = 1:10]

    if rank == root
        @test MPI.Reduce(send_arr, +, root, MPI.COMM_WORLD) ≈ [Double64(sz*i)/10 for i = 1:10] rtol=sz*eps(Double64)
    else
        @test MPI.Reduce(send_arr, +, root, MPI.COMM_WORLD) === nothing
    end

    MPI.Barrier( MPI.COMM_WORLD )
end

GC.gc()
MPI.Finalize()
@test MPI.Finalized()


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

mutable struct InstanceInfo <: SwaggerModel
    id::Union{ Nothing, String } # id
    state::Union{ Nothing, String } # state
    vmm_version::Union{ Nothing, String } # vmm_version

    function InstanceInfo(;id=nothing, state=nothing, vmm_version=nothing)
        o = new()
        set_field!(o, :id, id)
        set_field!(o, :state, state)
        set_field!(o, :vmm_version, vmm_version)
        o
    end
end # type InstanceInfo

const _name_map_InstanceInfo = Dict{String,Symbol}(["id"=>:id, "state"=>:state, "vmm_version"=>:vmm_version])
const _field_map_InstanceInfo = Dict{Symbol,String}([:id=>"id", :state=>"state", :vmm_version=>"vmm_version"])
Swagger.name_map(::Type{ InstanceInfo }) = _name_map_InstanceInfo
Swagger.field_map(::Type{ InstanceInfo }) = _field_map_InstanceInfo

const _allowed_InstanceInfo_state = ["Uninitialized", "Starting", "Running", "Halting", "Halted"]

function check_required(o::InstanceInfo)
    (o.id === nothing) && (return false)
    (o.vmm_version === nothing) && (return false)
    true
end

function validate_field(o::InstanceInfo, name::Symbol, val)
    if name === :state
        Swagger.validate_param(name, "InstanceInfo", :enum, val, _allowed_InstanceInfo_state)
    end
end


text:
nums = open("inputs/day09.txt") do f
    parse.(Int64, readlines(f))
end

function solve(nums, step=25)::Union{Int64,Nothing}
    start = 1
    for i in start + step:length(nums)
        slice = view(nums, i:(i + step))
        last = nums[i + step]
        search = function ()
            for a in slice
                for b in slice
                    if last == a + b
                        return true
                    end
                end
            end
            return false
        end
        if !search()
            return last
        end
    end
end

function solve2(nums, target)::Union{Int64,Nothing}
    for start in 1:length(nums)
        mx = nums[start]
        mn = nums[start]
        acc = nums[start]
        for i in (start + 1):length(nums)
            acc += nums[i]
            mx = max(mx, nums[i])
            mn = min(mn, nums[i])
            if acc > target
                break
            end
            if acc == target
                return mn + mx
            end
        end
    end
end

function main()
    part1 = solve(nums)
    println("Part 1: ", part1)
    println("Part 2: ", solve2(nums, part1))
end

main()


text:
### A Pluto.jl notebook ###
# v0.16.2

using Markdown
using InteractiveUtils

# ╔═╡ 122d494e-326b-11ec-0e29-551537ccbc53
begin
	import Pkg
	Pkg.activate(".")
end

# ╔═╡ 0350a012-0b7d-4f14-89e1-d01b860b2928
begin
	using Revise
	using PlotlyBase
	using PlutoUtils
	using TelecomUtils
end

# ╔═╡ 07fbb59d-cb6c-4066-878d-130ba1951ee1
p = generate_hex_lattice(1;M=10)

# ╔═╡ 3f56bf85-5d42-4d2e-8d68-b2999ac507f9
col = TelecomUtils.generate_colors(p)

# ╔═╡ 4cf7dfc7-b846-48b4-a902-7c7cb21ee6cd
scatter(p;marker = attr(
	color = TelecomUtils.generate_colors(p,4),
	colorscale = "Jet",
),mode = "markers") |> Plot

# ╔═╡ Cell order:
# ╠═122d494e-326b-11ec-0e29-551537ccbc53
# ╠═0350a012-0b7d-4f14-89e1-d01b860b2928
# ╠═07fbb59d-cb6c-4066-878d-130ba1951ee1
# ╠═3f56bf85-5d42-4d2e-8d68-b2999ac507f9
# ╠═4cf7dfc7-b846-48b4-a902-7c7cb21ee6cd


text:
passwordMatcher = r"(\d*)-(\d*) (.): (\w*)"

validPasswords = 0
for line in eachline()
    result = match(passwordMatcher, line)
    firstIndex = parse(Int32, result[1])
    secondIndex = parse(Int32, result[2])
    char = result[3]
    password = result[4]

    if xor(password[firstIndex] == char[1], password[secondIndex] == char[1])
        global validPasswords += 1
    end
end

println(validPasswords)


text:
for (T1, T2) in ((:T, :(float(T))), (:(Union{Missing,T}), :(Union{Missing,float(T)})))
  @eval begin

    # unweighted windowed function application that tapers

    function running(fun2::Function, data1::V, data2::V, windowspan::Int) where {T, V<:AbstractVector{$T1}}
        ndata   = min(length(data1), length(data2))
        nvals   = nrolled(ndata, windowspan)
        ntapers = ndata - nvals

        result = zeros($T2, ndata)

        result[1:ntapers] = tapers2(fun2, data1[1:ntapers], data2[1:ntapers])
        ntapers += 1
        result[ntapers:ndata] = rolling(fun2, data1, data2, windowspan)

        if isnan(result[1]) && ndata > 1
            result[1] = result[2]
        end

        return result
    end

    function running(fun2::Function, data1::V, data2::V, windowspan::Int, first::A) where {A, T, V<:AbstractVector{$T1}}
        ndata   = min(length(data1), length(data2))
        nvals   = nrolled(ndata, windowspan)
        ntapers = ndata - nvals

        result = zeros($T2, ndata)

        result[1:ntapers] = tapers2(fun2, data1[1:ntapers], data2[1:ntapers])
        ntapers += 1
        result[ntapers:ndata] = rolling(fun2, data1, data2, windowspan)

        result[1] = T(first)

        return result
    end

    # weighted windowed function application that tapers

    function running(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::F) where
                     {T, N<:Number, V<:AbstractVector{$T1}, F<:Vector{N}}
        ndata   = min(length(data1), length(data2))
        nvals   = nrolled(ndata, windowspan)
        ntapers = ndata - nvals

        result = zeros($T2, ndata)

        result[1:ntapers] = tapers2(fun2, data1[1:ntapers], data2[1:ntapers], weighting[end-(ntapers-1):end])
        ntapers += 1
        result[ntapers:ndata] = rolling(fun2, data1, data2, windowspan, weighting)

        if isnan(result[1]) && ndata > 1
            result[1] = result[2]
        end

        return result
    end

    function running(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::F, first::A) where
                     {A, T, V<:AbstractVector{$T1}, N<:Number, F<:Vector{N}}

        ndata   = min(length(data1), length(data2))
        nvals   = nrolled(ndata, windowspan)
        ntapers = ndata - nvals

        result = zeros($T2, ndata)

        result[1:ntapers] = tapers2(fun2, data1[1:ntapers], data2[1:ntapers], weighting[end-(ntapers-1):end])
        ntapers += 1
        result[ntapers:ndata] = rolling(fun2, data1, data2, windowspan, weighting)

        result[1] = T(first)

        return result
    end

    running(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::W) where
                    {T, V<:AbstractVector{$T1}, N<:Number, W<:AbstractWeights} =
        running(fun2, data1, data2, windowspan, weighting.values)

    running(fun2::Function, data1::V, data2::V, windowspan::Int, weighting::W, first::A) where
                    {T, V<:AbstractVector{$T1}, A, N<:Number, W<:AbstractWeights} =
        running(fun2, data1, data2, windowspan, weighting.values, first)

  end
end