text:
"""
    measurement_counts(samples::Matrix{Pair{String, Int}}; fillzeros::Bool = true)

Generate a dictionary containing the measurement counts for a set
of projectors, given a set of single-shot samples with different
measurement bases (i.e. QST).
"""
function measurement_counts(samples::Matrix{Pair{String, Int}}; fillzeros::Bool = true)
  counts = Dict{Tuple,Dict}()
  N = size(samples,2)

  if N > 8
    error("Full QST restricted to N ≤ 8")
  end
  # Loop over each measurement in the data
  for n in 1:size(samples,1)
    # Extract the measurement basis and  and outcome
    basis    = Tuple(first.(samples[n,:]))
    outcome  = Tuple(last.(samples[n,:]))

    # If new basis, add dictionary
    if !haskey(counts,basis)
      counts[basis] = Dict{Tuple,Int64}()
    end
    # Record outcome
    if !haskey(counts[basis],outcome)
      counts[basis][outcome] = 0
    end
    counts[basis][outcome] += 1
  end

  if fillzeros
    # Fill the counts with zeros for bitstring never observed
    for (k1,v1) in counts
      for index in 0:1<<N-1
        bitstring = reverse(Tuple(digits.(index,base=2,pad=N)'))
        if !haskey(counts[k1],bitstring)
          counts[k1][bitstring] = 0
        end
      end
    end
  end
  return counts
end

measurement_counts(samples::Matrix{String}; kwargs...) =
  measurement_counts(convertdatapoints(samples); kwargs...)

"""
    measurement_counts(data::Matrix{Pair{String,Pair{String, Int}}}; fillzeros::Bool = true)

Generate a dictionary containing the measurement counts for a set
of input states and measurement projectors (i.e. QPT).
"""
function measurement_counts(data::Matrix{Pair{String,Pair{String, Int}}}; fillzeros::Bool = true)
  if size(data,2) > 8
    error("Full QPT restricted to N ≤ 4")
  end
  newdata = []
  input_states = first.(data)
  measurements = last.(data)
  bases = first.(measurements)
  outcomes = last.(measurements)
  for n in 1:size(bases,1)
    m = convertdatapoint(outcomes[n,:], bases[n,:])
    tmp = []
    for j in 1:size(bases,2)
      push!(tmp, input_states[n,j])
      push!(tmp, m[j])
    end
    push!(newdata, tmp)
  end
  return measurement_counts(permutedims(hcat(newdata...)); fillzeros = fillzeros)
end

"""
    empirical_probabilities(counts::Dict{Tuple,Dict})

Compute the probabilities from a set of measurement counts.
"""
function empirical_probabilities(counts::Dict{Tuple,Dict})
  probs = Dict{Tuple,Dict{Tuple,Float64}}(deepcopy(counts))

  for basis in keys(probs)
    # Get count dictionary for 2^N projectors in this basis
    counts_in_basis = probs[basis]

    # Total number of counts
    tot_counts = sum(values(counts_in_basis))

    # Build probabilities
    for projector in keys(counts_in_basis)
      probs[basis][projector] /= tot_counts
    end
  end
  return probs
end

empirical_probabilities(samples::Array; fillzeros::Bool=true) =
  empirical_probabilities(measurement_counts(samples;fillzeros=fillzeros))

"""
    projector_matrix(probs::AbstractDict; process::Bool = false, return_probs::Bool = false)

Return the projector matrix, where each row corresponds to the vectorized
projector into different measurement bases contained into an input
probability dictionary.

If `return_probs=true`, return also the 1d vector of probabilities (i.e. the
1D flattening of the dictionary).
"""
function design_matrix(probs::AbstractDict; process::Bool = false, return_probs::Bool = false)
  A = []
  p̂ = []

  for (basis,projectors) in probs
    for (outcome,probability) in projectors
      Π_list = []
      for j in 1:length(outcome)
        if process && isodd(j)
          Πj = 0.5 * transpose(gate("Id") + (1-2*outcome[j]) * gate(basis[j]))
        else
          Πj = 0.5 * (gate("Id") + (1-2*outcome[j]) * gate(basis[j]))
        end
        push!(Π_list,Πj)
      end
      Π = reduce(kron,Π_list)
      push!(A,conj(vec(Π)))
      push!(p̂,probability)
    end
  end

  return_probs && return copy(transpose(hcat(A...))), float.(p̂)
  return copy(transpose(hcat(A...)))
end

"""
    make_PSD(ρ::AbstractArray{<:Number,2})

Transform a non-positive density operator into positive.
The output (positive) density matrix is obtained by minimizing
the 2-norm distance between the input density matrix and a positive
density matrix.
"""
function make_PSD(ρ::AbstractArray{<:Number,2})
  μ, μ⃗ = eigen(ρ; sortby=x->-real(x))
  a = 0.0
  μ = real(μ)
  λ = zeros(length(μ))
  i = length(μ)
  while (μ[i] + a/i < 0.0)
    a += μ[i]
    λ[i] = 0.0
    i -= 1
  end
  for j in reverse(1:i)
    λ[j] = μ[j] + a/i
  end

  ρ_PSD = zeros(ComplexF64,length(μ),length(μ))
  for k in 1:length(μ)
    ρ_PSD += λ[k] * (μ⃗[:,k] * μ⃗[:,k]')
  end
  return ρ_PSD
end


text:
function crop!(df, tStart::DateTime, tStop::DateTime)
  flt1 = df[:date] .>= tStart
  flt2 = df[:date] .<= tStop
  flt = flt1 & flt2
  df = df[flt, :]
end

function crop!(df, tStart::ASCIIString, tStop::ASCIIString)
  crop!(df, DateTime(tStart), DateTime(tStop))
end

function crop(df, tStart::DateTime, tStop::DateTime)
  flt1 = df[:date] .>= tStart
  flt2 = df[:date] .<= tStop
  flt = flt1 & flt2
  return df[flt, :]
end

function crop(df, tStart::ASCIIString, tStop::ASCIIString)
  crop(df, DateTime(tStart), DateTime(tStop))
end

function clean!(df)
  flt = ones(Bool, size(df, 1))
  for var in [:peakArea, :peakCenter, :peakWidth]
    flt &= get_outliers(df, var=var)
  end
  df = df[flt, :]
end

function clean(df)
  flt = ones(Bool, size(df, 1))
  for var in [:peakArea, :peakCenter, :peakWidth]
    flt &= get_outliers(df, var=var)
  end
  try
    flt_offNadir = df[:offNadirAngle_deg] .<= 25.0
    flt &= flt_offNadir
  end
  return df[flt, :]
end

function clean_h5(path)
  fileNames = readdir(path)
  for fname in fileNames
    fid = h5open(joinpath(path, fname),"r+")
    for grp in fid
      for dataset in grp
        if ismatch(r"M06.2", name(dataset))
          println(name(dataset))
          o_delete(grp, name(dataset))
        end
      end
    end
    close(fid)
  end
end

function doCO2Correction(tCO, nCO, tCO2, nCO2)
  dtMax = DateTime(2000,1,1,0,10) - DateTime(2000,1,1,0,0)
  nCOCorrected = Float64[]
  tCOCorrected = DateTime[]
  wasCorrected = zeros(Int64, length(nCO))
  for i=1:length(nCO)
    for j=1:length(nCO2)
      if ((tCO2[j] > tCO[i]) & ((tCO2[j] - tCO[i]) < dtMax))
        newValue = nCO[i] - 0.099 * nCO2[j]
        if newValue > 0
          push!(nCOCorrected, newValue)
          push!(tCOCorrected, tCO[i])
          wasCorrected[i] = 1
          break
        end
      end
    end
  end
  return tCOCorrected, nCOCorrected, wasCorrected
end

function moving_avg(y, window_size)
  myAvg = copy(y)
  for i=window_size+1:length(y)-window_size
    myAvg[i] = median(y[i-window_size:i+window_size])
  end
  return myAvg
end

function pds2h5(path, homeDir="/home/abieler/rosetta/data/dfms/h5/")
  fileNames = readdir(path)

  for fileName in sort(fileNames)
    # remove the .TAB of fileName
    bName = fileName[1:end-4]
    yrStr, tStr, modeStr = matchall(r"(\d+)", bName)
    if (modeStr[end] == '2') & (modeStr != "0572") & (modeStr[1:2] != "06")
      h5Str = joinpath(homeDir, yrStr[1:6] * ".h5")
      y, gainStep, t, m0, p_pt = parseDataFileBothRows(joinpath(path, fileName))
      mstr = "m" * string(round(Int,m0))
      dataSetName = mstr * "/" * bName
      try
        h5write(h5Str, joinpath(dataSetName, "gainStep"), gainStep)
        h5write(h5Str, joinpath(dataSetName, "commandedMass"), m0)
        h5write(h5Str, joinpath(dataSetName, "p_pt"), p_pt)
        h5write(h5Str, joinpath(dataSetName, "y"), y)
        h5write(h5Str, joinpath(dataSetName, "tStart"), string(t))
      catch
      end
    end
  end
end

function remove_outliers(df; N=50, var=:peakArea)
  flt = get_outliers(df, N=N, var=var)
  return df[flt,:]
end

function remove_outliers!(df; N=50, var=:peakArea)
  flt = get_outliers(df, N=N, var=var)
  df = df[flt,:]
end

function get_outliers(df; N=50, var=:peakArea)
  y = convert(Array, df[var])
  nDataPoints = length(y)
  flt = ones(Bool, nDataPoints)
  for i = N+1:nDataPoints-N
    data = append!(y[i-N:i-1], y[i+1:i+N])
    if abs(median(data) - y[i]) > (3*std(data))
      flt[i] = false
    end
  end
  return flt
end

ext: jl


text:
function generate_basic()
    # Waveguide definition
    segment_ranges = [0, 500e3, 1000e3, 1500e3]
    hprimes = [72.0, 74, 75, 76]
    betas = [0.28, 0.30, 0.32, 0.35]
    b_mags = fill(50e-6, length(segment_ranges))
    b_dips = fill(π/2, length(segment_ranges))
    b_azs = fill(0.0, length(segment_ranges))
    ground_sigmas = [0.001, 0.001, 0.0005, 0.0001]
    ground_epsrs = [4, 4, 10, 10]

    # Transmitter
    frequency = 24e3

    # Outputs
    output_ranges = collect(0:20e3:2000e3)

    input = ExponentialInput()
    input.name = "basic"
    input.description = "Test ExponentialInput"
    input.datetime = Dates.now()

    input.segment_ranges = segment_ranges
    input.hprimes = hprimes
    input.betas = betas
    input.b_mags= b_mags
    input.b_dips = b_dips
    input.b_azs = b_azs
    input.ground_sigmas = ground_sigmas
    input.ground_epsrs = ground_epsrs
    input.frequency = frequency
    input.output_ranges = output_ranges

    json_str = JSON3.write(input)

    open("basic.json","w") do f
        write(f, json_str)
    end

    return nothing
end

function generate_table()
    # Waveguide definition
    segment_ranges = [0, 500e3, 1000e3, 1500e3]

    hprimes = [72.0, 74, 75, 76]
    betas = [0.28, 0.30, 0.32, 0.35]

    altitude = 40e3:110e3
    density = [Vector{Float64}(undef, length(altitude)) for i = 1:length(segment_ranges)]
    collision_frequency = similar(density)
    for i in eachindex(segment_ranges)
        density[i] = waitprofile.(altitude, hprimes[i], betas[i])
        collision_frequency[i] = electroncollisionfrequency.(altitude)
    end

    b_mags = fill(50e-6, length(segment_ranges))
    b_dips = fill(π/2, length(segment_ranges))
    b_azs = fill(0.0, length(segment_ranges))
    ground_sigmas = [0.001, 0.001, 0.0005, 0.0001]
    ground_epsrs = [4, 4, 10, 10]

    # Transmitter
    frequency = 24e3

    # Outputs
    output_ranges = collect(0:20e3:2000e3)

    input = TableInput()
    input.name = "table"
    input.description = "Test TableInput"
    input.datetime = Dates.now()

    input.segment_ranges = segment_ranges
    input.altitude = altitude
    input.density = density
    input.collision_frequency = collision_frequency
    input.b_mags= b_mags
    input.b_dips = b_dips
    input.b_azs = b_azs
    input.ground_sigmas = ground_sigmas
    input.ground_epsrs = ground_epsrs
    input.frequency = frequency
    input.output_ranges = output_ranges

    json_str = JSON3.write(input)

    open("table.json","w") do f
        write(f, json_str)
    end

    return nothing
end

function generate_batchtable()
    N = 2
    rep(v) = repeat(v, 1, N)

    # Waveguide definition
    nsegments = 4
    segment_ranges = [0, 500e3, 1000e3, 1500e3]
    b_mags = fill(50e-6, nsegments)
    b_dips = fill(π/2, nsegments)
    b_azs = fill(0.0, nsegments)
    ground_sigmas = [0.001, 0.001, 0.0005, 0.0001]
    ground_epsrs = [4, 4, 10, 10]

    hprimes = [72.0, 74, 75, 76]
    betas = [0.28, 0.30, 0.32, 0.35]

    altitude = 40e3:110e3
    density = [Vector{Float64}(undef, length(altitude)) for i = 1:length(segment_ranges)]
    collision_frequency = similar(density)
    for i in eachindex(segment_ranges)
        density[i] = waitprofile.(altitude, hprimes[i], betas[i])
        collision_frequency[i] = electroncollisionfrequency.(altitude)
    end

    # Transmitter
    frequency = 24e3

    # Outputs
    output_ranges = collect(0:20e3:2000e3)

    binput = BatchInput{TableInput}()
    binput.name = "batchtable"
    binput.description = "Test BatchInput with TableInput"
    binput.datetime = Dates.now()

    inputs = Vector{TableInput}(undef, N)
    for i in eachindex(inputs)
        input = TableInput()
        input.name = "$i"
        input.description = "Test TableInput $i"
        input.datetime = Dates.now()

        input.segment_ranges = segment_ranges
        input.altitude = altitude
        input.density = density
        input.collision_frequency = collision_frequency
        input.b_mags= b_mags
        input.b_dips = b_dips
        input.b_azs = b_azs
        input.ground_sigmas = ground_sigmas
        input.ground_epsrs = ground_epsrs
        input.frequency = frequency
        input.output_ranges = output_ranges

        inputs[i] = input
    end

    binput.inputs = inputs

    json_str = JSON3.write(binput)

    open("batchtable.json","w") do f
        write(f, json_str)
    end

    return nothing
end

function generate_batchbasic()
    N = 2
    rep(v) = repeat(v, 1, N)

    # Waveguide definition
    nsegments = 4
    segment_ranges = [0, 500e3, 1000e3, 1500e3]
    b_mags = fill(50e-6, nsegments)
    b_dips = fill(π/2, nsegments)
    b_azs = fill(0.0, nsegments)
    ground_sigmas = [0.001, 0.001, 0.0005, 0.0001]
    ground_epsrs = [4, 4, 10, 10]

    hprimes = rep([72.0, 74, 75, 76])
    betas = rep([0.28, 0.30, 0.32, 0.35])

    # Transmitter
    frequency = 24e3

    # Outputs
    output_ranges = collect(0:20e3:2000e3)

    binput = BatchInput{ExponentialInput}()
    binput.name = "batchbasic"
    binput.description = "Test BatchInput with ExponentialInput"
    binput.datetime = Dates.now()

    inputs = Vector{ExponentialInput}(undef, N)
    for i in eachindex(inputs)
        input = ExponentialInput()

        input.name = "$i"
        input.description = "ExponentialInput $i"
        input.datetime = binput.datetime
        input.segment_ranges = segment_ranges
        input.hprimes = hprimes[:,i]
        input.betas = betas[:,i]
        input.b_mags = b_mags
        input.b_dips = b_dips
        input.b_azs = b_azs
        input.ground_sigmas = ground_sigmas
        input.ground_epsrs = ground_epsrs
        input.frequency = frequency
        input.output_ranges = output_ranges

        inputs[i] = input
    end

    binput.inputs = inputs

    json_str = JSON3.write(binput)

    open("batchbasic.json","w") do f
        write(f, json_str)
    end

    return nothing
end

function read_basic()
    open("basic.json","r") do f
        v = JSON3.read(f, ExponentialInput)
        return v
    end
end

function read_table()
    open("table.json", "r") do f
        v = JSON3.read(f, TableInput)
        return v
    end
end

function read_batchbasic()
    open("batchbasic.json","r") do f
        v = JSON3.read(f, BatchInput{ExponentialInput})
        return v
    end
end

function read_corrupted_basic()
    # missing the ground_sigmas field
    open("corrupted_basic.json","r") do f
        v = JSON3.read(f, ExponentialInput)
        return v
    end
end

function run_basic()
    s = LMP.parse("basic.json")
    output = LMP.buildrun(s)
    return output
end

function run_table()
    s = LMP.parse("table.json")
    output = LMP.buildrun(s)
    return output
end

function run_batchbasic()
    s = LMP.parse("batchbasic.json")
    output = LMP.buildrun(s)
    return output
end

function run_batchbasicsave()
    s = LMP.parse("batchbasic.json")
    output = LMP.buildrunsave("batchbasictest.json", s)
    sres = LMP.parse("batchbasictest.json")
    return sres
end

function run_batchtablesave()
    s = LMP.parse("batchtable.json")
    output = LMP.buildrunsave("batchtabletest.json", s)
    sres = LMP.parse("batchtabletest.json")
    return sres
end

function basic_lmp()
    output = LMP.propagate("basic.json")

    # Compare to a `propagate` call with waveguide arguments
    tx = Transmitter(24e3)
    rx = GroundSampler(0:20e3:2000e3, Fields.Ez)

    # From `generate_basic()`
    segment_ranges = [0, 500e3, 1000e3, 1500e3]
    hprimes = [72.0, 74, 75, 76]
    betas = [0.28, 0.30, 0.32, 0.35]
    b_mags = fill(50e-6, length(segment_ranges))
    b_dips = fill(π/2, length(segment_ranges))
    b_azs = fill(0.0, length(segment_ranges))
    ground_sigmas = [0.001, 0.001, 0.0005, 0.0001]
    ground_epsrs = [4, 4, 10, 10]

    wvg = SegmentedWaveguide(
        [HomogeneousWaveguide(
            BField(b_mags[i], b_dips[i], b_azs[i]),
            Species(QE, ME, z->waitprofile(z, hprimes[i], betas[i]; cutoff_low=40e3), electroncollisionfrequency),
            Ground(ground_epsrs[i], ground_sigmas[i]),
            segment_ranges[i]
        ) for i = 1:4]
    )

    _, a, p = propagate(wvg, tx, rx)

    @test isapprox(a, output.amplitude, atol=0.1)
    @test isapprox(p, output.phase, atol=deg2rad(1))

    return output
end

@testset "IO.jl" begin
    @info "Testing IO"

    generate_basic()
    generate_table()
    generate_batchbasic()
    generate_batchtable()

    @test LMP.iscomplete(read_basic())
    @test LMP.validlengths(read_basic())
    @test LMP.iscomplete(read_corrupted_basic()) == false
    @test LMP.parse("basic.json") isa ExponentialInput
    @test LMP.parse("basic.json", ExponentialInput) isa ExponentialInput

    @test LMP.iscomplete(read_table())
    @test LMP.validlengths(read_table())
    @test LMP.parse("table.json") isa TableInput
    @test LMP.parse("table.json", TableInput) isa TableInput

    @test LMP.iscomplete(read_batchbasic())
    @test LMP.validlengths(read_batchbasic())
    @test LMP.parse("batchbasic.json") isa BatchInput{ExponentialInput}

    @info "  Running:"
    @info "    Segmented Wait ionospheres..."
    @test run_basic() isa BasicOutput
    @info "    Segmented tabular ionospheres..."
    @test run_table() isa BasicOutput
    @info "    Multiple segmented Wait ionospheres..."
    @test run_batchbasic() isa BatchOutput{BasicOutput}
    @info "    Multiple segmented Wait ionospheres, appending..."
    @test run_batchbasicsave() isa BatchOutput{BasicOutput}
    @info "    Multiple segmented tabular ionospheres, appending."
    @test run_batchtablesave() isa BatchOutput{BasicOutput}

    @info "    `propagate` segmented Wait ionospheres"
    @test basic_lmp() isa BasicOutput

    isfile("basic.json") && rm("basic.json")
    isfile("table.json") && rm("table.json")
    isfile("batchbasic.json") && rm("batchbasic.json")
    isfile("batchtable.json") && rm("batchtable.json")
    isfile("batchbasictest.json") && rm("batchbasictest.json")
    isfile("batchtabletest.json") && rm("batchtabletest.json")
    isfile("basic_output.json") && rm("basic_output.json")
end


text:
#Parent type of all simulations
abstract type @par(AbstractSimulation) end
#=
Simulation will encapsule different structs for:
  Time-stepping, PassiveScalar, DensityStratification, LESModel, Forcing Scheme
=#

#Traits
@par haspassivescalar(s::Type{T}) where {T<:@par(AbstractSimulation)} =
    PassiveScalarType !== NoPassiveScalar
haspassivescalar(s::AbstractSimulation) = haspassivescalar(typeof(s))

@par hasdensity(s::Type{T}) where {T<:@par(AbstractSimulation)}  =
    DensityStratificationType !== NoDensityStratification
hasdensity(s::AbstractSimulation) = hasdensity(typeof(s))

@par hasles(s::Type{T}) where {T<:@par(AbstractSimulation)}  =
    LESModelType !== NoLESModel
hasles(s::AbstractSimulation) = hasles(typeof(s))

@par hasdensityles(s::Type{T}) where {T<:@par(AbstractSimulation)}  =
    (hasdensity(T) && hasdensityles(DensityStratificationType))
hasdensityles(s::AbstractSimulation) = hasdensityles(typeof(s))

@par haspassivescalarles(s::Type{T}) where {T<:@par(AbstractSimulation)}  =
    (haspassivescalar(T) && haspassivescalarles(PassiveScalarType))
haspassivescalarles(s::AbstractSimulation) = haspassivescalarles(typeof(s))

@par hasforcing(s::Type{T}) where {T<:@par(AbstractSimulation)}  =
    ForcingType !== NoForcing
hasforcing(s::AbstractSimulation) = hasforcing(typeof(s))

@par hashyperviscosity(s::Type{T}) where {T<:@par(AbstractSimulation)} =
    HyperViscosityType !== NoHyperViscosity
hashyperviscosity(s::AbstractSimulation) = hashyperviscosity(typeof(s))

@par is_vorticityEquation(s::Type{T}) where {T<:@par(AbstractSimulation)} =
    EquationType === VorticityEquation
is_vorticityEquation(s::AbstractSimulation) = is_vorticityEquation(typeof(s))

@par need_vorticity(s::Type{T}) where {T<:@par(AbstractSimulation)} =
    is_vorticityEquation(T) | hasles(T)
need_vorticity(s::AbstractSimulation) = need_vorticity(typeof(s))

struct @par(Simulation) <: @par(AbstractSimulation)
    u::VectorField{Float64,3,2,false}
    rhs::VectorField{Float64,3,2,false}
    reductionh::Vector{Float64}
    reductionv::Vector{Float64}
    reductionh_fil::Vector{Float64}
    reductionv_fil::Vector{Float64}
    xspec::PaddedArray{Float64,1,0,true}
    yspec::PaddedArray{Float64,1,0,true}
    equation::EquationType
    timestep::VelocityTimeStepType
    passivescalar::PassiveScalarType
    densitystratification::DensityStratificationType
    lesmodel::LESModelType
    forcing::ForcingType
    hyperviscosity::HyperViscosityType
    iteration::Base.RefValue{Int}
    time::Base.RefValue{Float64}
    dtoutput::Int
    dtstats::Int
    dtspec::Int
    hspec::Array{Float64,3}
    vspec::Array{Float64,3}
    specH::Array{Float64,1}
    tspecH::Array{Float64,1}
    specV::Array{Float64,1}
    kspec::ScalarField{Float64,3,2,false}
    nlstats::Base.RefValue{Tuple{Float64,Float64,Float64,Float64}}
    pressstats::Base.RefValue{Tuple{Float64,Float64,Float64,Float64}}

    @par function @par(Simulation)(u::VectorField,equation,timestep,passivescalar,densitystratification,lesmodel,forcing,hv,iteration,time,dtout,dtstats,dtspecs)

        rhs = similar(u)
        kspec = ScalarField{Float64}(size(u.r),(LX,LY,LZ))

        reductionh = zeros(THR ? Threads.nthreads() : 1)
        reductionv = zeros(THR ? Threads.nthreads() : 1)

        reductionh_fil = !STAT_FIL ? zeros(1) : zeros(THR ? Threads.nthreads() : 1)
        reductionv_fil = !STAT_FIL ? zeros(1) : zeros(THR ? Threads.nthreads() : 1)

        xspec = PaddedArray(NRX)
        yspec = PaddedArray(NY)

        hspec = zeros(Float64,size(u))
        vspec = zeros(Float64,size(u))
        specH = zeros(Float64,length(KH))
        tspecH = zeros(Float64,length(KH))
        specV = zeros(Float64,length(KRZ))
        nlstats = Ref((1.0,1.0,1.0,1.0))
        pressstats = Ref((1.0,1.0,1.0,1.0))

        return @par(new)(u,rhs,reductionh,reductionv,reductionh_fil,reductionv_fil,xspec,yspec,equation,timestep,passivescalar,densitystratification,lesmodel,forcing,hv,iteration,time,dtout,dtstats,dtspecs,hspec,vspec,specH,tspecH,specV,kspec,nlstats,pressstats)
    end

end

is_output_time(s) = (s.timestep.x.iteration[] != 0) & (mod(s.iteration[],s.dtoutput) == 0)
is_stats_time(s) = (mod(s.iteration[],s.dtstats) == 0)
is_spec_time(s) = (mod(s.iteration[],s.dtspec) == 0)

##################################### Equation type #########################################

struct VorticityEquation end

struct NavierStokesEquation{T}
    uu::SymTrTenField{T,3,2,false}
end

msg(::VorticityEquation) = "Navier-Stokes in rotational form ((∇×u)×u)"
msg(::NavierStokesEquation) = "Navier-Stokes in conservation form (∇⋅(uu))"

##################################### Equation type #########################################

@inline @par nuh(s::Type{T}) where {T<:@par(AbstractSimulation)} = nuh(HyperViscosityType)
@inline nuh(s::AbstractSimulation) = nuh(typeof(s))

@inline @par get_hyperviscosity_exponent(s::Type{T}) where {T<:@par(AbstractSimulation)} = get_hyperviscosity_exponent(HyperViscosityType)
@inline get_hyperviscosity_exponent(s::AbstractSimulation) = get_hyperviscosity_exponent(typeof(s))

@par function Base.show(io::IO,s::@par(Simulation))
smsg = """
Fluid Flow Simulation
$(msg(s.equation))

nx: $(NRX)
ny: $NY
nz: $NZ
x domain size: $(LX)*2π
y domain size: $(LY)*2π
z domain size: $(LZ)*2π

Kinematic Viscosity: $(ν)

Velocity time-stepping method: $(typeof(s.timestep.x))
Dealias cutoff: $(DEALIAS_TYPE)
Number of Threads: $NT
"""
smsg = join((smsg,msg.(getfield.(Ref(s),sim_fields))...))#msg(s.passivescalar),
#  msg(s.densitystratification),
#  msg(s.lesmodel),
#  msg(s.forcing)))

print(io,smsg)
end

#LES models
include("LESmodels/les_types.jl")

# Simulaiton with Scalar fields ===================================================================================================================================================

abstract type AbstractPassiveScalar{L,TT,α,dρdz,Gdirec} end

    haspassivescalarles(a::Type{T}) where {L,TT,α,dρdz,Gdirec,T<:AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}} =
        L === NoLESScalar ? false : true
    @inline haspassivescalarles(a::AbstractPassiveScalar) = haspassivescalarles(typeof(a))

    has_les_scalar_vorticity_model(a::Type{T}) where {L,TT,α,dρdz,Gdirec,T<:AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}} =
        is_vorticity_model(L)
    @inline has_les_scalar_vorticity_model(a::AbstractPassiveScalar) = has_les_scalar_vorticity_model(typeof(a))
    @inline @par has_les_scalar_vorticity_model(s::Type{T}) where {T<:@par(AbstractSimulation)} = has_les_scalar_vorticity_model(PassiveScalarType)

    diffusivity(a::Type{T}) where {L,TT,α,dρdz,Gdirec,T<:AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}} =
        α
    @inline diffusivity(a::AbstractPassiveScalar) = diffusivity(typeof(a))

    meangradient(a::Type{T}) where {L,TT,α,dρdz,Gdirec,T<:AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}} =
        dρdz
    @inline meangradient(a::AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}) where {L,TT,α,dρdz,Gdirec} = meangradient(typeof(a))

    graddir(a::Type{T}) where {L,TT,α,dρdz,Gdirec,T<:AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}} =
        Gdirec
    @inline graddir(a::AbstractPassiveScalar{L,TT,α,dρdz,Gdirec}) where {L,TT,α,dρdz,Gdirec} = graddir(typeof(a))

    initialize!(a::AbstractPassiveScalar,s::AbstractSimulation) = initialize!(a.timestep,parent(real(a.rhs)),diffusivity(a),s)

    statsheader(a::AbstractPassiveScalar) = "scalar,scalarp2,dscalardxp2,dscalardyp2,dscalardzp2,scalardiss"

    stats(a::AbstractPassiveScalar,s::AbstractSimulation) = scalar_stats(a.φ,a,s)

struct NoPassiveScalar <: AbstractPassiveScalar{NoLESScalar,nothing,nothing,nothing,nothing} end

    initialize!(a::NoPassiveScalar,s::AbstractSimulation) = nothing

    statsheader(a::NoPassiveScalar) = ""

    stats(a::NoPassiveScalar,s::AbstractSimulation) = (),()

    msg(a::NoPassiveScalar) = "\nPassive Scalar: No passive scalar\n"

struct PassiveScalar{L,TTimeStep, α #=Difusitivity = ν/Pr =#,
                  dρdz #=Linear mean profile=#, Gdirec #=Axis of mean profile =#} <: AbstractPassiveScalar{L,TTimeStep,α,dρdz,Gdirec}
    φ::ScalarField{Float64,3,2,false}
    rhs::ScalarField{Float64,3,2,false}
    flux::VectorField{Float64,3,2,false}
    lesmodel::L
    timestep::TTimeStep
    reduction::Vector{Float64}

    function PassiveScalar{TT,α,dρdz,Gdirec}(ρ,les,timestep) where {TT,α,dρdz,Gdirec}
        ρrhs = similar(ρ)
        flux = VectorField(size(ρ.field.r),(ρ.kx.l,ρ.ky.l,ρ.kz.l))
        reduction = zeros(THR ? Threads.nthreads() : 1)
        return new{typeof(les),TT,α,dρdz,Gdirec}(ρ,ρrhs,flux,les,timestep,reduction)
    end
end

msg(a::PassiveScalar{L,TT,α,dρdz,Gdirec}) where {L,TT,α,dρdz,Gdirec} = """

Passive Scalar: true
Scalar diffusivity: $(α)
Scalar mean gradient: $(dρdz)
Scalar mean gradient direction: $(Gdirec)
Scalar time-stepping method: $(TT)
Scalar LES model: $(L)

"""
# ==========================================================================================================
# Hyper viscosity Type

include("Hyperviscosity/hyperviscosity_types.jl")

# ==========================================================================================================
# Density Stratification Types

include("density_stratification_types.jl")

# ==========================================================================================
# Forcing Scheme

include("Forcing_methods/forcing_types.jl")

# Initializan function =========================================================================================================================================================================================================

function parameters(d::Dict)

    nx = parse(Int,d[:nx])
    ny = parse(Int,d[:ny])
    nz = parse(Int,d[:nz])

    ncx = div(nx,2)+1
    lcs = ncx*ny*nz
    lcv = 3*lcs
    lrs = 2*lcs
    lrv = 2*lcv

    lx = Float64(eval(Meta.parse(d[:xDomainSize])))
    ly = Float64(eval(Meta.parse(d[:yDomainSize])))
    lz = Float64(eval(Meta.parse(d[:zDomainSize])))
    ν = Float64(eval(Meta.parse(d[:kinematicViscosity])))

    start = haskey(d,:start) ? parse(Int,d[:start]) : 0
    starttime = haskey(d,:startTime) ? parse(Float64,d[:startTime]) : 0.0
    dtstat = parse(Int,d[:dtStat])
    dtspecs = parse(Int,d[:dtSpec])
    dtout = parse(Int,d[:writeTime])

    if haskey(d,:equation)
        if d[:equation] == "NS"
            equation = NavierStokesEquation{Float64}(SymTrTenField((nx,ny,nz),(lx,ly,lz)))
        else
            equation = VorticityEquation()
        end
    else
        equation = VorticityEquation()
    end

    @info("Reading initial velocity field u1.$start u2.$start u3.$start")
    u = VectorField{Float64}("u1.$start","u2.$start","u3.$start")

    kxp = reshape(rfftfreq(nx,lx),(ncx,1,1))
    kyp = reshape(fftfreq(ny,ly),(1,ny,1))
    kzp = reshape(fftfreq(nz,lz),(1,1,nz))

    haskey(d,:threaded) ? (tr = parse(Bool,d[:threaded])) : (tr = true)

    tr || FFTW.set_num_threads(1)
    nt = tr ? Threads.nthreads() : 1

    b = Globals.splitrange(lrs, nt)

    haskey(d,:cutoff) ? (cutoffr = Float64(eval(Meta.parse(d[:cutoff])))) : (cutoffr = 2/3)

    cutoff = (cutoffr*kxp[end])^2

    if haskey(d,:hyperViscosity)
        if d[:hyperViscosity] == "spectralBarrier"
            initkp = haskey(d,:initk) ? parse(Float64,d[:initk]) : cutoffr*2KX[end] / 3
            endK = haskey(d,:endk) ? parse(Float64,d[:endk]) : sqrt(cutoffr*3)*KX[end]
            hyperviscositytype = SpectralBarrier(initkp,endK)
        else
            νh = parse(Float64,d[:hyperViscosityCoefficient])
            m = haskey(d,:hyperViscosityM) ? parse(Int,d[:hyperViscosityM]) : 2
            hyperviscositytype = HyperViscosity{νh,m}()
        end
    else
        hyperviscositytype = NoHyperViscosity()
    end

    integrator = haskey(d,:velocityTimeStep) ? Symbol(d[:velocityTimeStep]) : :ETD3rdO
    variableTimeStep = haskey(d,:variableTimestepFlag) ? parse(Bool,d[:variableTimestepFlag]) : true
    cfl = haskey(d,:cfl) ? parse(Float64,d[:cfl]) : 0.18
    idt = haskey(d,:dt) ? Float64(eval(Meta.parse(d[:dt]))) : 0.0
    vtimestep = if integrator === :Euller
        VectorTimeStep{cfl}(Euller{variableTimeStep,idt}(Ref(idt)),Euller{variableTimeStep,idt}(Ref(idt)),Euller{variableTimeStep,idt}(Ref(idt)))
    elseif integrator === :Adams_Bashforth3rdO
        VectorTimeStep{cfl}(Adams_Bashforth3rdO{variableTimeStep,idt}(),Adams_Bashforth3rdO{variableTimeStep,idt}(),Adams_Bashforth3rdO{variableTimeStep,idt}())
    else
        VectorTimeStep{cfl}(ETD3rdO(variableTimeStep,idt,hyperviscositytype,ν),ETD3rdO(variableTimeStep,idt,hyperviscositytype,ν),ETD3rdO(variableTimeStep,idt,hyperviscositytype,ν))
    end

    if haskey(d,:passiveScalar)
        α = ν/Float64(eval(Meta.parse(d[:scalarPr])))
        dρdz = Float64(eval(Meta.parse(d[:scalarGradient])))
        @info("Reading initial scalar field scalar.$start")
        rho = isfile("scalar.$start") ? ScalarField{Float64}("scalar.$start") : ScalarField{Float64}((nx,ny,nz),(lx,ly,lz))
        scalardir = haskey(d,:scalarDirection) ? Symbol(d[:scalarDirection]) : :z
        scalartimestep = if integrator === :Euller
            Euller{variableTimeStep,idt}(Ref(idt))
        elseif integrator === :Adams_Bashforth3rdO
            Adams_Bashforth3rdO{variableTimeStep,idt}()
        else
            ETD3rdO(variableTimeStep,idt,hyperviscositytype,α)
        end

        lestypescalar = if haskey(d,:lesModel)
            if haskey(d,:scalarLesModel)
                if d[:scalarLesModel] == "vorticityDiffusion"
                    c = haskey(d,:scalarVorticityCoefficient) ? parse(Float64,d[:scalarVorticityCoefficient]) : 1.0
                    VorticityDiffusion(c)
                else
                    EddyDiffusion()
                end
            else
                EddyDiffusion()
            end
        else
            NoLESScalar()
        end

        scalartype = PassiveScalar{typeof(scalartimestep),α,dρdz,scalardir}(rho,lestypescalar,scalartimestep)
    else
        scalartype = NoPassiveScalar()
    end

    densitytype = stratificationtype(d,start,integrator,nx,ny,nz,lx,ly,lz,variableTimeStep,idt)

    lestype = les_types(d,nx,ny,nz,lx,ly,lz)

    forcingtype = forcing_model(d,nx,ny,nz,ncx,lx)

    s = Simulation{typeof(equation),typeof(vtimestep),
        typeof(scalartype),typeof(densitytype),typeof(lestype),typeof(forcingtype),
        typeof(hyperviscositytype)}(u,equation,vtimestep,scalartype,densitytype,lestype,forcingtype,hyperviscositytype,Ref(start),Ref(starttime),dtout,dtstat,dtspecs)
  #
    return s
end

parameters() = parameters(readglobal())


text:
using Mongo

type TwMongoDB <: TwitterDB
  db::MongoClient
  col::MongoCollection

  function TwMongoDB(db_name::AbstractString, col_name::AbstractString)
    mongo_client = MongoClient()
    mongo_collection = MongoCollection(mongo_client, db_name, col_name)
    new(mongo_client, mongo_collection)
  end

end

function add_entry(tmc::TwMongoDB, d::Dict{AbstractString, Any})
  insert(tmc.col, d)
end

function close_db(tmc::TwMongoDB)
  #Not needed
end

max_stars_repo_name: nick-paul/JuliaTwitter

max_issues_repo_name: nick-paul/JuliaTwitter

max_forks_repo_name: nick-paul/JuliaTwitter


text:
using Krylov
using LinearAlgebra
using SparseArrays
using CUDA

function square_preconditioned(n :: Int=10)
        A   = ones(n, n) + (n-1) * Matrix(I, n, n)
        b   = 10.0 * [1:n;]
        M⁻¹ = 1/n * Matrix(I, n, n)
        return A, b, M⁻¹
end

A, b, N = square_preconditioned(100)
spA = sparse(A)
spN = sparse(N)
y, status = Krylov.bicgstab(spA, copy(b), N=spN)
@show length(status.residuals)
@show y

cuspA = CUDA.CUSPARSE.CuSparseMatrixCSR(spA)
cuspN = CUDA.CUSPARSE.CuSparseMatrixCSR(spN)
cub = CuVector(b)
cuy, status = Krylov.bicgstab(cuspA, cub, N=cuspN)
@show cuy
@show length(status.residuals)


text:
"""
Definition of simulation data types for juliaPBF
    All variables are treated using the struct
    in SimulationDataStruct. the result of phy
    sical calculation will be updated into the
    struct.
"""
## EXPORT TYPES
export
    ParticleData,
    BoundaryData,
    GridData
## EXPORT TYPE FUNCTIONs
export
    # Change
    changeGridID,
    changeLambda,
    changeVel,
    changePos,
    changeTemppos,
    changeMass,
    # Add
    addVel,
    addPos,
    addTemppos,
    addmMass,
    # ETC
    updatePos

""" a Particle Data Type """
struct ParticleData ## Array of Struct
        gridID      ::Int64
        vel         ::Vec2
        pos         ::Vec2
        temppos     ::Vec2
        lambda      ::Float64
        mass        ::Float64
        phase       ::Int64 # SimulationData와 맞춰진. 이걸타고 들어가서 Fluid/Solid 나뉠것임
end
# CHANGE VALUEs
changeGridID(a::ParticleData, b::Int64)   = ParticleData(b,        a.vel,   a.pos,   a.temppos,   a.lambda,   a.mass,   a.phase)
changeLambda(a::ParticleData, b::Float64) = ParticleData(a.gridID, a.vel,   a.pos,   a.temppos,   b,          a.mass,   a.phase)
changeVel(a::ParticleData, b::Vec2)       = ParticleData(a.gridID, b,       a.pos,   a.temppos,   a.lambda,   a.mass,   a.phase)
changePos(a::ParticleData, b::Vec2)       = ParticleData(a.gridID, a.vel,   b,       a.temppos,   a.lambda,   a.mass,   a.phase)
changeTemppos(a::ParticleData, b::Vec2)   = ParticleData(a.gridID, a.vel,   a.pos,   b,           a.lambda,   a.mass,   a.phase)
changeMass(a::ParticleData, b::Float64)   = ParticleData(a.gridID, a.vel,   a.pos,   a.temppos,   a.lambda,   b,        a.phase)
# ADD VALUEs
addVel(a::ParticleData, b::Vec2)          = ParticleData(a.gridID, a.vel+b, a.pos,   a.temppos,   a.lambda,   a.mass,   a.phase)
addPos(a::ParticleData, b::Vec2)          = ParticleData(a.gridID, a.vel,   a.pos+b, a.temppos,   a.lambda,   a.mass,   a.phase)
addTemppos(a::ParticleData, b::Vec2)      = ParticleData(a.gridID, a.vel,   a.pos,   a.temppos+b, a.lambda,   a.mass,   a.phase)
addMass(a::ParticleData, b::Float64)      = ParticleData(a.gridID, a.vel,   a.pos,   a.temppos,   a.lambda,   a.mass+b, a.phase)
# ETC
updatePos(a::ParticleData)                = ParticleData(a.gridID, a.vel,   a.temppos,a.temppos,  a.lambda,   a.mass,   a.phase)
isless(a::ParticleData,b::ParticleData)   = isless(a.gridID::Int64,b.gridID::Int64) #for Sorting by gridID

""" Boundary Data Type """
struct BoundaryData
    points          ::Vector{Vec2}
    lines           ::Vector{iVec2} # Point Index 짝궁
end
""" Grid Data Type """
struct GridData
    nNearGrid       ::Vector{Int64} # Cumulate 형식
    nearGridIndex   ::Vector{Int64} # Near Grid Index를 1줄행렬로 가지고 있음
    nGridParticles  ::Vector{Int64} # 그리드 내 속한 입자 수
end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Group operation action" begin
    G = GroupManifold(NotImplementedManifold(), Manifolds.MultiplicationOperation())
    A_left = GroupOperationAction(G)
    A_right = GroupOperationAction(G, RightAction())

    types = [Matrix{Float64}]

    @test group_manifold(A_left) === G
    @test base_group(A_left) == G
    @test repr(A_left) == "GroupOperationAction($(repr(G)), LeftAction())"
    @test repr(A_right) == "GroupOperationAction($(repr(G)), RightAction())"

    @test switch_direction(LeftAction()) == RightAction()
    @test switch_direction(RightAction()) == LeftAction()

    for type in types
        a_pts = convert.(type, [reshape(i:(i + 3), 2, 2) for i in 1:3])
        m_pts = convert.(type, [reshape((i + 2):(i + 5), 2, 2) for i in 1:3])

        atol = eltype(m_pts[1]) == Float32 ? 1e-5 : 1e-10

        test_action(
            A_left,
            a_pts,
            m_pts;
            test_optimal_alignment=false,
            test_diff=false,
            atol=atol,
        )

        test_action(
            A_right,
            a_pts,
            m_pts;
            test_optimal_alignment=false,
            test_diff=false,
            atol=atol,
        )
    end

    G = SpecialOrthogonal(3)
    M = Rotations(3)
    A_left = GroupOperationAction(G)
    A_right = GroupOperationAction(G, RightAction())

    p = Matrix{Float64}(I, 3, 3)
    aω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
    mω = [2 .* ω for ω in aω]
    a_pts = [exp(M, p, hat(M, p, ωi)) for ωi in aω]
    m_pts = [exp(M, p, hat(M, p, ωi)) for ωi in mω]
    X_pts = [
        hat(M, p, [-1.0, 2.0, 0.5]),
        hat(M, p, [2.0, 1.0, 0.5]),
        hat(M, p, [0.5, 0.5, 0.5]),
    ]

    @test group_manifold(A_left) === G
    @test base_group(A_left) == G
    @test repr(A_left) == "GroupOperationAction($(repr(G)), LeftAction())"
    @test repr(A_right) == "GroupOperationAction($(repr(G)), RightAction())"

    test_action(A_left, a_pts, m_pts, X_pts; test_optimal_alignment=true, test_diff=true)

    test_action(A_right, a_pts, m_pts, X_pts; test_optimal_alignment=true, test_diff=true)
end


text:
# Base.:+(a::Nothing,b::Nothing) = nothing
# Base.:-(a::Nothing,b::Nothing) = nothing

function simple_allowed_reactions(
    dgs::Vector{<:Union{Missing,Real}}, #this shouldn't contain `nothing`s though
    threshold::Real,
    allow_missing::Bool=false)

    allowed_reactions = dgs .<= threshold
    if allow_missing==true
        allowed_reactions[ismissing.(dgs)] .= true
    else
        allowed_reactions[ismissing.(dgs)] .= false
    end

    Vector{Bool}(allowed_reactions)
end

"""
Indicates which reactions can proceed based on the dg threshold.

Returns two Vector{Bool} of length equal to rids, one for forward reactions and one for backwards.

Note: `allow_unbalanced==false` will only allow reactions with `is_balanced==nothing`
      when `allow_nothings==true`.
      `allow_unbalanced==true` will only allow reactions with `is_balanced==nothing`
      when `allow_nothings==true`
"""
function filter_reactions_by_dg(
    threshold::Real,
    env_key::Any,
    rids::IDs,
    rstructs::Reactions;
    allow_nothings::Bool=false,
    allow_unbalanced::Bool=false,
    allow_within_ci::Bool=true)

    mean_dgs = Union{Nothing,Missing,Real}[rstructs[rid].metadata["dg"][env_key]["standard_dg_prime_value"] for rid in rids]
    cis = Union{Nothing,Missing,Real}[rstructs[rid].metadata["dg"][env_key]["standard_dg_prime_ci"] for rid in rids]
    ## convert to missing because it can handle operators
    mean_dgs[isnothing.(mean_dgs)] .= missing
    cis[isnothing.(cis)] .= missing

    ## identify best case scenario dg values
    if allow_within_ci==true
        forward_dgs = Vector{Union{Missing,Real}}(mean_dgs - cis)
        backward_dgs = Vector{Union{Missing,Real}}(-mean_dgs - cis)
    else
        forward_dgs = Vector{Union{Missing,Real}}(mean_dgs)
        backward_dgs = Vector{Union{Missing,Real}}(-mean_dgs)
    end

    forward_allowed = simple_allowed_reactions(forward_dgs,threshold,allow_nothings)
    backward_allowed = simple_allowed_reactions(backward_dgs,threshold,allow_nothings)

    ## gather balance data
    balances = Union{Nothing,Missing,Bool}[rstructs[rid].metadata["dg"][env_key]["is_balanced"] for rid in rids]
    ## convert to missing because it can handle operators
    balances[isnothing.(balances)] .= missing

    if (allow_unbalanced == true) & (allow_nothings == true)
        ## don't use balance data
        return forward_allowed, backward_allowed

    elseif (allow_unbalanced == true) & (allow_nothings == false)
        forward_allowed[ismissing.(balances)] .= false
        backward_allowed[ismissing.(balances)] .= false
        return forward_allowed, backward_allowed

    elseif (allow_unbalanced == false) & (allow_nothings == true)
        balances[ismissing.(balances)] .= true

        forward_allowed[balances .== false] .= false
        backward_allowed[balances .== false] .= false
        return forward_allowed, backward_allowed

    elseif (allow_unbalanced == false) & (allow_nothings == false)
        balances[ismissing.(balances)] .= false

        forward_allowed[balances .== false] .= false
        backward_allowed[balances .== false] .= false
        return forward_allowed, backward_allowed
    end
end

max_stars_repo_name: hbsmith/BioXP

max_issues_repo_name: hbsmith/BioXP

max_forks_repo_name: hbsmith/BioXP


text:
module Falcons

using ReferenceFrameRotations
using Healpix
using LinearAlgebra
using Base.Threads
using StaticArrays
using ProgressMeter

include("./function/func4scan.jl")
include("./function/scanning.jl")
include("./function/mapmake.jl")
include("./function/make_scanmap.jl")
#include("./function/focalplane.jl")

export ScanningStrategy, gen_ScanningStrategy, get_pointings, get_pointing_pixels, period2rpm, get_pointings_xyz_tuple, get_pointings_tuple
export TwoTelescopes_ScanningStrategy2map, TwoTelescopes_ScanningStrategy2MapInfo, ThreeTelescopes_ScanningStrategy2map
export ScanningStrategy2map, Mapmaking, pixtod2hitmap, angtod2hitmap, array2map, xyztod2hitmap
export get_psiDataBase, ScanningStrategy2MapInfo
#export lft_focalplane_configration, mft_focalplane_configration, hft_focalplane_configration, get_FP_each_freq
#export pickup_wafer, pickup_freq, ang2xy, longitude2pix, ang2xy_projected_hft, ang2xy_projected_mft

end


text:
@timedtestset "addition" begin
    @timedtestset "Correlated GPs" begin
        rng, N, N′, D, gpc = MersenneTwister(123456), 5, 6, 2, GPC()
        X, X′ = ColVecs(randn(rng, D, N)), ColVecs(randn(rng, D, N′))
        f1 = wrap(GP(1, SEKernel()), gpc)
        f2 = wrap(GP(2, SEKernel()), gpc)
        f3 = f1 + f2
        f4 = f1 + f3
        f5 = f3 + f4

        for (n, (fp, fa, fb)) in enumerate([(f3, f1, f2), (f4, f1, f3), (f5, f3, f4)])
            Σp = cov(fa(X)) + cov(fb(X)) + cov(fa(X), fb(X)) + cov(fb(X), fa(X))
            ΣpXX′ = cov(fa(X), fa(X′)) + cov(fb(X), fb(X′)) + cov(fa(X), fb(X′)) +
                cov(fa(X), fb(X′))
            @test mean(fp(X)) ≈ mean(fa(X)) + mean(fb(X))
            @test cov(fp(X)) ≈ Σp
            @test cov(fp(X), fp(X′)) ≈ ΣpXX′
            @test cov(fp(X′), fp(X)) ≈ transpose(ΣpXX′)
            @test cov(fp(X), fa(X′)) ≈ cov(fa(X), fa(X′)) + cov(fb(X), fa(X′))
            @test cov(fp(X′), fa(X)) ≈ cov(fa(X′), fa(X)) + cov(fb(X′), fa(X))
            @test cov(fa(X), fp(X′)) ≈ cov(fb(X), fa(X′)) + cov(fa(X), fa(X′))
            @test cov(fa(X′), fp(X)) ≈ cov(fb(X′), fa(X)) + cov(fa(X′), fa(X))
        end

        @timedtestset "Consistency Tests" begin
            P, Q = 4, 3
            x0, x1, x2, x3 = randn(rng, P), randn(rng, Q), randn(rng, Q), randn(rng, P)
            abstractgp_interface_tests(f3, f1, x0, x1, x2, x3)
            abstractgp_interface_tests(f2 - f1, f1, x0, x1, x2, x3)
        end
    end
    @timedtestset "Verify mean / kernel numerically" begin
        rng, N, D = MersenneTwister(123456), 5, 6
        X = ColVecs(randn(rng, D, N))
        c, f = randn(rng), wrap(GP(5, SEKernel()), GPC())

        @test mean((f + c)(X)) == mean(f(X)) .+ c
        @test mean((f + c)(X)) == c .+ mean(f(X))
        @test cov((f + c)(X)) == cov(f(X))
        @test cov((c + f)(X)) == cov(f(X))

        @test mean((f - c)(X)) == mean(f(X)) .- c
        @test mean((c - f)(X)) == c .- mean(f(X))
        @test cov((f - c)(X)) == cov(f(X))
        @test cov((c - f)(X)) == cov(f(X))

        x = randn(rng, N + D)
        @test mean((f + sin)(x)) == mean(f(x)) + map(sin, x)
        @test mean((sin + f)(x)) == map(sin, x) + mean(f(x))
        @test cov((f + sin)(x)) == cov(f(x))
        @test cov((sin + f)(x)) == cov(f(x))

        @timedtestset "Consistency Tests" begin
            P, Q = 5, 3
            x0, x1, x2, x3 = randn(rng, P), randn(rng, Q), randn(rng, Q), randn(rng, P)
            abstractgp_interface_tests(c + f, f, x0, x1, x2, x3)
        end
    end
    @timedtestset "Standardised Tests (independent sum)" begin
        standard_1D_tests(
            MersenneTwister(123456),
            Dict(:l1=>0.5, :l2=>2.3),
            θ->begin
                gpc = GPC()
                f1 = θ[:l1] * wrap(GP(sin, SEKernel()), gpc)
                f2 = θ[:l2] * wrap(GP(cos, SEKernel()), gpc)
                f3 = f1 + f2
                return f3, f3
            end,
            collect(range(-1.5, 1.5; length=5)),
            collect(range(-1.0, 0.5; length=3)),
        )
    end
    @timedtestset "Standardised Tests (correlated sum)" begin
        standard_1D_tests(
            MersenneTwister(123456),
            Dict(:l1=>0.5, :l2=>2.3),
            θ->begin
                gpc = GPC()
                f1 = θ[:l1] * wrap(GP(sin, SEKernel()), gpc)
                f2 = θ[:l2] * wrap(GP(cos, SEKernel()), gpc)
                f3 = f1 + f2
                f4 = f1 + f3
                f5 = f3 + f4
                return f5, f5
            end,
            collect(range(-1.5, 1.5; length=5)),
            collect(range(-1.0, 0.5; length=3)),
        )
    end
end

# # θ = Dict(:l1=>0.5, :l2=>2.3);
# x, z = collect(range(-5.0, 5.0; length=512)), collect(range(-5.0, 5.0; length=128));
# y = rand(GP(sin, SqExponentialKernel(), GPC())(x, 0.1));

# foo_logpdf = (x, y) -> begin
#     gpc = GPC()
#     f = GP(sin, SqExponentialKernel(), gpc)
#     return logpdf(f(x, 0.1), y)
# end

# foo_elbo = (x, y, z) -> begin
#     f = GP(0, SqExponentialKernel(), GPC())
#     return elbo(f(x, 0.1), y, f(z, 0.001))
# end

# @benchmark foo_logpdf($x,  $y)
# @benchmark Zygote.pullback(foo_logpdf, $x, $y)

# let
#     z, back = Zygote.pullback(foo_logpdf, x, y)
#     @benchmark $back($(randn()))
# end

# let
#     foo = function(x, y)
#         fx = GP(0, SqExponentialKernel(), GPC())(x, 0.1)
#         C = cholesky(Symmetric(cov(fx)))
#         return logdet(C) + Xt_invA_X(C, y)
#     end
#     display(@benchmark Zygote.pullback($foo, $x, $y))
#     z_pw, back_pw = Zygote.pullback(foo, x, y)
#     @benchmark $back_pw(randn())
# end

# @benchmark foo_elbo($x, $y, $z)
# @benchmark Zygote.pullback(foo_elbo, $x, $y, $z)

# let
#     L, back = Zygote.pullback(foo_elbo, x, y, z)
#     @benchmark $back($L)
# end

# θ->begin
#     gpc = GPC()
#     f1 = GP(sin, SqExponentialKernel(l=θ[:l1]), gpc)
#     f2 = GP(cos, SqExponentialKernel(l=θ[:l2]), gpc)
#     f3 = f1 + f2
#     return f3, f3
# end,
# 13, 11,


text:
using CausalityToolsBase
using DelayEmbeddings

"""
    te_embed(source::AbstractArray{T, 1}, target::AbstractArray{T, 1},
        k = 1, l = 1, m = 1, n = 0, η = 1, τ = 1)

Create an appropriate custom delay reconstruction from `source` and `target` for the
purpose of computing transfer entropy from `source` to `target`.

# Delay reconstruction

Denote the source time series by `s(t)` and target time series by `t(t)`. Then the
embedding vectors are on the the form

``(T_{f}, T_{pp}, S_{pp})`` where

``T_f = (t(t + k*η), ..., t(t + 2*η), t(t + η))``,
``T_pp = (t(t), t(t - τ), t(t - 2*τ), ..., t(t - lτ))``,
``S_pp = (s(t), s(t - τ), s(t - 2*τ), ..., s(t - mτ))``

and the total dimension of the reconstructed space is `k + l + m`.

# Arguments

- **`source`**: A source data series.
- **`target`**: A target data series.
- **`k`**: The dimension of the ``T_f`` component of the delay reconstruction.
- **`l`**: The dimension of the ``T_{pp}`` component of the delay reconstruction.
- **`m`**: The dimension of the ``S_{pp}`` component of the delay reconstruction.

# Keyword arguments
- **`η`**: The prediction lag that goes into the ``T_f`` component of the
    delay reconstruction.
- **`τ`**: The embedding lag that goes into the ``T_{pp}`` and ``S_{pp}``
    components of the delay reconstruction.

# Examples

```julia
x, y, z = rand(100), rand(100), rand(100)

# Get embedding points and instructions for computing marginals.

# TE(x -> y)
# ==================
pts, vars = te_embed(x, y) # default k = l = m = 1
pts, vars = te_embed(x, y, 1, 1, 1) # specify k, l and m

pts, vars = te_embed(x, y, η = 2) # specify prediction lag
pts, vars = te_embed(x, y, τ = 2) # specify embedding lag

# Compute TE
# ======================
transferentropy(pts, vars, RectangularBinning(5), VisitationFrequency())
```
"""
function te_embed(source::AbstractArray{T, 1}, target::AbstractArray{T, 1},
        k = 1, l = 1, m = 1; η = 1, τ = 1) where T

    length(source) == length(target) || throw(ArgumentError("length(source) != length(target)"))
    k >= 1 || throw(ArgumentError("k=$k must be >= 1 (so that Tf has entries)"))
    l >= 1 || throw(ArgumentError("l=$l must be >= 1 (so that Tpp has entries)"))
    m >= 1 || throw(ArgumentError("m=$m must be >= 1 (so that Spp has entries)"))
    k + l + m >= 3 || throw(ArgumentError("k + l + m must be at least 3 for meaningful TE computations"))

    # Positions (ks and ls both refer to the target time series,
    # while ms refer to the source time series).
    pos_ks = repeat([2], k)
    pos_ls = repeat([2], l)
    pos_ms = repeat([1], m)
    pos = Positions([pos_ks; pos_ls; pos_ms])

    # Lags
    lags_ks = [k*η for k in k:-1:1]
    lags_ls = [-l*τ for l in 0:(l - 1)]
    lags_ms = [-m*τ for m in 0:(m - 1)]

    lags = Lags([lags_ks; lags_ls; lags_ms])

    v = TEVars(Tf = 1:k |> collect, Tpp = 1+(k):l+(k) |> collect, Spp = 1+(l+k):m+(l+k) |> collect)
    customembed(Dataset(source, target), pos, lags), v
end

"""
    te_embed(source, target, cond, k = 1, l = 1, m = 1, n = 1; η = 1, τ = 1)

Create an appropriate custom delay reconstruction from `source`, `target` and `cond` for the
purpose of computing transfer entropy from `source` to `target` conditioned on `cond`.

# Delay reconstruction

Denote the entries in the source time series by `s(t)`, the entries of the target time
series by `t(t)` and the entries of the conditional time series by `c(t)`. Then the embedding
vectors are on the the form

``(T_{f}, T_{pp}, S_{pp}, C_{pp})`` where

``T_f = (t(t + k*η), ..., t(t + 2*η), t(t + η))``,
``T_pp = (t(t), t(t - τ), t(t - 2*τ), ..., t(t - lτ))``,
``S_pp = (s(t), s(t - τ), s(t - 2*τ), ..., s(t - mτ))``,
``C_pp = (c(t), c(t - τ), c(t - 2*τ), ..., c(t - mτ))``

and the total dimension of the reconstructed space is `k + l + m + n`.

# Arguments

- **`source`**: A source data series.
- **`target`**: A target data series.
- **`cond`**: A data series on which to condition.
- **`k`**: The dimension of the ``T_f`` component of the delay reconstruction.
- **`l`**: The dimension of the ``T_{pp}`` component of the delay reconstruction.
- **`m`**: The dimension of the ``S_{pp}`` component of the delay reconstruction.
- **`n`**: The dimension of the ``C_{pp}`` component of the delay reconstruction.

# Keyword arguments
- **`η`**: The prediction lag that goes into the ``T_f`` component of the
    delay reconstruction.
- **`τ`**: The embedding lag that goes into the ``T_{pp}``, ``S_{pp}`` and ``C_{pp}``
    components of the delay reconstruction.

# Examples

```julia
x, y, z = rand(100), rand(100), rand(100)

# Get embedding points and instructions for computing marginals.

# Conditional TE(x -> y | z)
# ==========================
pts, vars = te_embed(x, y, z) # default k = l = m = n = 1
pts, vars = te_embed(x, y, z, η = 2, τ = 3) # specify prediction and embedding lag
pts, vars = te_embed(x, y, 2, 1, 2, 1) # specify k = 2, l = 1, m = 2, n = 1

# Compute TE
# ======================
transferentropy(pts, vars, RectangularBinning(5), VisitationFrequency())
```
"""
function te_embed(source::AbstractArray{T, 1},
            target::AbstractArray{T, 1},
            cond::AbstractArray{T, 1},
            k = 1, l = 1, m = 1, n = 1; η = 1, τ = 1) where T

    length(source) == length(target) || throw(ArgumentError("length(source) != length(target)"))
    k >= 1 || throw(ArgumentError("k=$k must be >= 1 (so that Tf has entries)"))
    l >= 1 || throw(ArgumentError("l=$l must be >= 1 (so that Tpp has entries)"))
    m >= 1 || throw(ArgumentError("m=$m must be >= 1 (so that Spp has entries)"))
    n >= 1 || throw(ArgumentError("n=$n must be >= 1 (so that Cpp has entries)"))

    k + l + m + n >= 4 || throw(ArgumentError("k + l + m + n must be at least 4 for meaningful TE computations"))

    # Positions (ks and ls both refer to the target time series,
    # while ms refer to the source time series).
    pos_ks = repeat([2], k)
    pos_ls = repeat([2], l)
    pos_ms = repeat([1], m)
    pos_ns = repeat([3], n)

    pos = Positions([pos_ks; pos_ls; pos_ms; pos_ns])

    # Lags
    lags_ks = [k*η for k in k:-1:1]
    lags_ls = [-l*τ for l in 0:(l - 1)]
    lags_ms = [-m*τ for m in 0:(m - 1)]
    lags_ns = [-n*τ for n in 0:(n - 1)]

    lags = Lags([lags_ks; lags_ls; lags_ms; lags_ns])

    v = TEVars(
        Tf = 1:k |> collect,
        Tpp = 1+(k):l+(k) |> collect,
        Spp = 1+(l+k):m+(l+k) |> collect,
        Cpp = 1+(l+k+m):n+(l+k+m) |> collect)

    customembed(Dataset(source, target, cond), pos, lags), v
end

export te_embed


text:
for I in sectorlist
    println("------------------------------------")
    println("Sector $I")
    println("------------------------------------")
    ti = time()
    @testset "Sector $I: Basic properties" begin
        s = (randsector(I), randsector(I), randsector(I))
        @test eval(Meta.parse(sprint(show,s[1]))) == s[1]
        @test @constinferred(hash(s[1])) == hash(deepcopy(s[1]))
        @test @constinferred(one(s[1])) == @constinferred(one(I))
        @constinferred dual(s[1])
        @constinferred dim(s[1])
        @constinferred frobeniusschur(s[1])
        @constinferred Nsymbol(s...)
        @constinferred Rsymbol(s...)
        @constinferred Bsymbol(s...)
        @constinferred Fsymbol(s..., s...)
        it = @constinferred s[1] ⊗ s[2]
        @constinferred ⊗(s..., s...)
        for i = 1:3
            @test 1 == @constinferred Nsymbol(s[i], conj(s[i]), one(s[i]))
        end
    end
    @timedtestset "Sector $I: Value iterator" begin
        @test eltype(values(I)) == I
        sprev = one(I)
        for (i, s) in enumerate(values(I))
            @test !isless(s, sprev) # confirm compatibility with sort order
            if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
                @test_throws ArgumentError values(I)[i]
                @test_throws ArgumentError TensorKit.findindex(values(I), s)
            elseif hasmethod(Base.getindex, Tuple{typeof(values(I)), Int})
                @test s == @constinferred (values(I)[i])
                @test TensorKit.findindex(values(I), s) == i
            end
            sprev = s
            i >= 10 && break
        end
        @test one(I) == first(values(I))
        if Base.IteratorSize(values(I)) == Base.IsInfinite() && I <: ProductSector
            @test_throws ArgumentError TensorKit.findindex(values(I), one(I))
        elseif hasmethod(Base.getindex, Tuple{typeof(values(I)), Int})
            @test (@constinferred TensorKit.findindex(values(I), one(I))) == 1
            for s in smallset(I)
                @test (@constinferred values(I)[TensorKit.findindex(values(I), s)]) == s
            end
        end
    end
    @timedtestset "Sector $I: fusion tensor and F-move and R-move" begin
        for a in smallset(I), b in smallset(I)
            for c in ⊗(a,b)
                X1 = permutedims(fusiontensor(a, b, c), (2, 1, 3, 4))
                X2 = fusiontensor(b, a, c)
                l =dim(a)*dim(b)*dim(c)
                R = LinearAlgebra.transpose(Rsymbol(a, b, c))
                sz = (l, convert(Int, Nsymbol(a, b, c)))
                @test reshape(X1, sz) ≈ reshape(X2, sz) * R
            end
        end
        for a in smallset(I), b in smallset(I), c in smallset(I)
            for e in ⊗(a, b), f in ⊗(b, c)
                for d in intersect(⊗(e, c), ⊗(a, f))
                    X1 = fusiontensor(a, b, e)
                    X2 = fusiontensor(e, c, d)
                    Y1 = fusiontensor(b, c, f)
                    Y2 = fusiontensor(a, f, d)
                    @tensor f1[-1,-2,-3,-4] := conj(Y2[a,f,d,-4])*conj(Y1[b,c,f,-3])*
                                                X1[a,b,e,-1] * X2[e,c,d,-2]
                    f2 = Fsymbol(a,b,c,d,e,f)*dim(d)
                    @test isapprox(f1, f2; atol = 1000*eps(), rtol = 1000*eps())
                end
            end
        end
    end
    @timedtestset "Sector $I: Unitarity of F-move" begin
        for a in smallset(I), b in smallset(I), c in smallset(I)
            for d in ⊗(a,b,c)
                es = collect(intersect(⊗(a,b), map(dual, ⊗(c,dual(d)))))
                fs = collect(intersect(⊗(b,c), map(dual, ⊗(dual(d),a))))
                Fblocks = Vector{Any}()
                for e in es
                    for f in fs
                        Fs = Fsymbol(a,b,c,d,e,f)
                        push!(Fblocks, reshape(Fs, (size(Fs, 1)*size(Fs, 2), size(Fs, 3)*size(Fs, 4))))
                    end
                end
                F = hvcat(length(fs), Fblocks...)
                @test F'*F ≈ one(F)
            end
        end
    end
    @timedtestset "Sector $I: Pentagon equation" begin
        for a in smallset(I), b in smallset(I), c in smallset(I), d in smallset(I)
            for f in ⊗(a,b), h in ⊗(c,d)
                for g in ⊗(f,c), i in ⊗(b,h)
                    for e in intersect(⊗(g,d), ⊗(a,i))
                        @tensor p1[λ,μ,ν,κ,ρ,σ] := Fsymbol(f,c,d,e,g,h)[λ,μ,ν,τ]*
                                                    Fsymbol(a,b,h,e,f,i)[κ,τ,ρ,σ]
                        p2 = zero(p1)
                        for j in ⊗(b,c)
                            @tensor p2[λ,μ,ν,κ,ρ,σ]  += Fsymbol(a,b,c,g,f,j)[κ,λ,α,β]*
                                                        Fsymbol(a,j,d,e,g,i)[β,μ,τ,σ]*
                                                        Fsymbol(b,c,d,i,j,h)[α,τ,ν,ρ]
                        end
                        @test isapprox(p1, p2; atol = 1000*eps(), rtol = 1000*eps())
                    end
                end
            end
        end
    end
    @timedtestset "Sector $I: Hexagon equation" begin
        for a in smallset(I), b in smallset(I), c in smallset(I)
            for e in ⊗(c,a), g in ⊗(c,b)
                for d in intersect(⊗(e,b), ⊗(a,g))
                    @tensor p1[α,β,μ,ν] := Rsymbol(c,a,e)[α,λ]*
                                            Fsymbol(a,c,b,d,e,g)[λ,β,γ,ν]*
                                            Rsymbol(b,c,g)[γ,μ]
                    p2 = zero(p1)
                    for f in ⊗(a,b)
                        @tensor p2[α,β,μ,ν] += Fsymbol(c,a,b,d,e,f)[α,β,δ,σ]*
                                                Rsymbol(c,f,d)[σ,ψ]*
                                                Fsymbol(a,b,c,d,f,g)[δ,ψ,μ,ν]
                    end
                    @test isapprox(p1, p2; atol = 1000*eps(), rtol = 1000*eps())
                end
            end
        end
    end
    tf = time()
    printstyled("Finished sector $I tests in ",
                string(round(tf-ti; sigdigits=3)),
                " seconds."; bold = true, color = Base.info_color())
    println()
end


text:
"Check if a teerm is ground (no variables)"
is_ground(t::Term) = error("Not implemented.")
is_ground(c::Constant) = true
is_ground(v::Variable) = false
is_ground(l::Structure) = all(is_ground(e) for e in l.arguments)
is_ground(l::Literal) = all(is_ground(a) for a in l.arguments)
is_ground(n::Negation) = is_ground(n.literal)
is_ground(c::Conj) = all(is_ground(e) for e in c.elements)
is_ground(cl::Clause) = is_ground(cl.head) && is_ground(c.body)

"Get all vars in a term"
get_vars(t::Term) = error("Not implemented.")
get_vars(t::Constant) = Set{Variable}()
get_vars(t::Variable) = Set{Variable}([t])
get_vars(s::Structure) = union((get_vars(t) for t in s.arguments)...)
get_vars(l::Literal) = union((get_vars(t) for t in l.arguments)...)
get_vars(n::Negation) = get_vars(n.literal)
get_vars(c::Conj) = union((get_vars(e) for e in c.elements)...)
get_vars(cl::Clause) = union(get_vars(cl.head), get_vars(cl.body))

"Get arity"
get_arity(c::Constant) = 1
get_arity(f::Functor) = f.arity
get_arity(p::Predicate) = p.arity
get_arity(s::Structure) = s.functor.arity
get_arity(l::Literal) = l.predicate.arity
get_arity(p::LPair) = 2
get_arity(v::Variable) = error("Cannot determine the arity of a variable")
get_arity(cl::Clause) = error("Cannot determine the arity of a clause")

"get functor name"
get_functor_name(f::Functor) = f.name
get_functor_name(p::Predicate) = p.name
get_functor_name(s::Structure) = s.functor.name
get_functor_name(l::Literal) = l.predicate.name
get_functor_name(c::Constant) = c.name
get_functor_name(v::Variable) = error("Variables do not have functors")
get_functor_name(e::Union{Clause,Conj}) = error("Not implemented")

"substitute variables"
subs(b::Bool, sub::Dict{Variable,Variable}) = b
subs(c::Constant, sub::Dict{Variable,Variable}) = c
subs(v::Variable, sub::Dict{Variable,Variable}) = get(sub, v, v)
subs(s::Structure, sub::Dict{Variable, Variable}) = Structure(s.functor, [subs(t, sub) for t in s.arguments])
subs(p::Proposition, subs::Dict{Variable,Variable}) = p
subs(l::List, sub::Dict{Variable,Variable}) = List([subs(e, sub) for e in l.elements])
subs(p::LPair, sub::Dict{Variable,Variable}) = LPair(subs(p.head, sub), subs(p.tail, sub))
subs(l::Literal, sub::Dict{Variable,Variable}) = Literal(l.predicate, [subs(a, sub) for a in l.arguments])
subs(c::Conj, sub::Dict{Variable, Variable}) = Conj([subs(e, sub) for e in c.elements])
subs(cl::Clause, sub::Dict{Variable,Variable}) = Clause(subs(cl.head, sub), subs(cl.body, sub))


text:
function compensated_summation(x::T, y::T, e::T) where {T}
    local err::T
    local res::T

    err = e   + x
    res = err + y
    e   = err + (y - res)

    return (res, e)
end

function compensated_summation(x::Vector{T}, y::Vector{T}, e::Vector{T}) where {T}
    local err::Vector{T}
    local res::Vector{T}
    local tmp::Vector{T}

    err = e   .+ x
    res = err .+ y
    tmp = y   .- res
    e  .= err .+ tmp

    return (res, e)
end

function L2norm(x)
    local l2::eltype(x) = 0
    for xᵢ in x
        l2 += xᵢ^2
    end
    l2
end

function l2norm(x)
    sqrt(L2norm(x))
end


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "LALO21"
cases = ["W=0.01", "W=0.05", "W=0.1"]
SUITE[model] = BenchmarkGroup()

include("laubloomis.jl")
validation = []
final_width = []

# ----------------------------------------
#  Case 1: smaller initial states
# ----------------------------------------
W = 0.01
prob = laubloomis(W=W);
alg = TMJets(abstol=1e-11, orderT=7, orderQ=1, adaptive=true);

# warm-up run
sol_1 = solve(prob, T=T, alg=alg);
sol_1z = overapproximate(sol_1, Zonotope);

# verify that specification holds
property = ρ(e4, sol_1z) < 4.5
push!(validation, Int(property))

# width of final box
width = ρ(e4, sol_1z[end]) + ρ(-e4, sol_1z[end])
push!(final_width, trunc(width, digits=4))
println("width of final box, case $(cases[1]) : $width")

# benchmark
SUITE[model][cases[1]] = @benchmarkable solve($prob, T=$T, alg=$alg)

# ----------------------------------------
# Case 2: intermediate initial states
# ----------------------------------------
W = 0.05
prob = laubloomis(W=W)
alg = TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=false)

# warm-up run
sol_2 = solve(prob, T=T, alg=alg)
sol_2z = overapproximate(sol_2, Zonotope)

# verify that specification holds
property = ρ(e4, sol_2z) < 4.5
push!(validation, property ? 1 : 0)

# width of final box
width = ρ(e4, sol_2z[end]) + ρ(-e4, sol_2z[end])
push!(final_width, trunc(width, digits=3))
println("width of final box, case $(cases[2]): $width")

# benchmark
SUITE[model][cases[2]] = @benchmarkable solve($prob, T=$T, alg=$alg)

# ----------------------------------------
# Case 3: larger initial states
# ----------------------------------------
W = 0.1
prob = laubloomis(W=W)
alg = TMJets(abstol=1e-12, orderT=7, orderQ=1, adaptive=false)

# warm-up run
sol_3 = solve(prob, T=T, alg=alg)
sol_3z = overapproximate(sol_3, Zonotope)

# verify that specification holds
property = ρ(e4, sol_3z) < 5.0
push!(validation, Int(property))

# width of final box
width = ρ(e4, sol_3z[end]) + ρ(-e4, sol_3z[end])
push!(final_width, trunc(width, digits=3))
println("width of final box, case $(cases[3]) : $width")

# benchmark
SUITE[model][cases[3]] = @benchmarkable solve($prob, T=$T, alg=$alg)

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
    t = median(results[model][c]).time * 1e-9
    runtimes[c] = t
end

for (i, c) in enumerate(cases)
    print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c]), $(final_width[i])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

fig = Plots.plot()

Plots.plot!(fig, sol_3z, vars=(0, 4), linecolor="green", color=:green, alpha=0.8)
Plots.plot!(fig, sol_2z, vars=(0, 4), linecolor="blue",  color=:blue, alpha=0.8)
Plots.plot!(fig, sol_1z, vars=(0, 4), linecolor="yellow", color=:yellow, alpha=0.8,
    tickfont=font(30, "Times"), guidefontsize=45,
    xlab=L"t", #\raisebox{-0.5mm}{\textcolor{white}{.}}",
    ylab=L"x_4", # \raisebox{2mm}{\textcolor{white}{.}}",
    xtick=[0., 5., 10., 15., 20.], ytick=[1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5.],
    xlims=(0., 20.), ylims=(1.5, 5.02),
    bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
    size=(1000, 1000))

Plots.plot!(fig, x->x, x->4.5, 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)
Plots.plot!(fig, x->x, x->5., 0., 20., line=2, color="red", linestyle=:dash, legend=nothing)

savefig("ARCH-COMP21-JuliaReach-LaubLoomis.png")


text:
JC69() = JC69rel()
JC69(λ::F; safe::Bool=true) where F <: Float64 = JC69abs(λ, safe=safe)

function JC69(θ_vec::A;
              safe::Bool=true) where A <: AbstractArray
  if length(θ_vec) == 0
    return JC69rel()
  elseif length(θ_vec) == 1
    return JC69abs(θ_vec[1], safe=safe)
  else
    error("Parameter vector length incompatiable with absolute or relative rate form of substitution model")
  end
end

function JC69rel(θ_vec::A;
                 safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 0
      error("Incorrect parameter vector length")
    end
  end
  return JC69rel()
end

function JC69abs(θ_vec::A;
                 safe::Bool=true) where A <: AbstractArray
  if safe
    if length(θ_vec) != 1
      error("Incorrect parameter vector length")
    end
  end
  return  JC69abs(θ_vec[1], safe=safe)
end


text:
using CategoricalArrays
using DataFrames
using Documenter
using SimpleANOVA
using Test

@testset "SimpleANOVA Tests" begin
    DocMeta.setdocmeta!(
        SimpleANOVA,
        :DocTestSetup,
        :(using SimpleANOVA);
        recursive=true
    )

    # Ensures that jldoctests also pass.
    doctest(SimpleANOVA)

    include("test_anova.jl")
    include("test_pretests.jl")
    include("test_contrasts.jl")
    include("test_show.jl")
end


text:
using DandelionWebSockets.Proxy
using DandelionWebSockets.Proxy: stopproxy
using Sockets: TCPSocket
using MbedTLS

mutable struct WebSocketConnectionImpl <: WebSocketConnection
    # `logic_proxy` forwards commands to the `ClientProtocol` object, in its own coroutine.
    logic_proxy::Union{ClientProtocolProxy, Nothing}
    # `reader` reads frames from the server.
    reader::Union{ServerReader, Nothing}
    # `rng` is used for random generation of masks when sending frames.
    rng::AbstractRNG
    # `ponger` keeps track of when a pong response is expected from the server.
    ponger::AbstractPonger
    # `pinger` requests that the logic send ping frames to the server at regular intervals.
    pinger::AbstractPinger

    # Requirement
    # @5_3-2-2 Masking uses strong source of entropy
    #
    # This is now handled by `RandomDevice`.

    WebSocketConnectionImpl() = new(nothing,
                                    nothing,
                                    RandomDevice(),
                                    Ponger(3.0, misses=3),
                                    Pinger(5.0))
end

"""
A WebSocket client, used to connect to a server, and send messages.

Note: The keyword arguments in the constructor are primarily for testing.
"""
mutable struct WSClient <: AbstractWSClient
    # WebSocketConnection maintains the state for a single connection.
    connection::Union{WebSocketConnection, Nothing}
    # `do_handshake` is a function that performs a HTTP Upgrade to a WebSocket connection.
    handshake::WebSocketHandshake

    WSClient(; handshake::WebSocketHandshake = HTTPHandshake(RandomDevice(), WebSocketUpgradeAdapter())) =
        new(nothing, handshake)
end

"""
    tcpnodelay(::IO)

Sets the TCP_NODELAY flag on a socket. This is a separate function only for testing purposes. It
can be implemented with a more specific type if the flag makes no sense for another `IO` subtype.
"""
tcpnodelay(io::TCPSocket) = ccall(:uv_tcp_nodelay, Cint, (Ptr{Nothing}, Cint), io, 1)
tcpnodelay(io::MbedTLS.SSLContext) = ccall(:uv_tcp_nodelay, Cint, (Ptr{Nothing}, Cint), io.bio, 1)

maybewrapio(io::TCPSocket) :: IO = io
maybewrapio(io::MbedTLS.SSLContext) :: IO = TLSBufferedIO(io)

"""
    connection_result(::WSClient, ::GoodHandshake, ::WebSocketHandler, fix_small_message_latency::Bool)

For a valid handshake, start all background tasks for this connection. This includes tasks for
reading and writing from the socket, as well as a task for user callback.
"""
function connection_result_(client::WSClient,
                            result::GoodHandshake,
                            handler::WebSocketHandler,
                            fix_small_message_latency::Bool)
    # Requirement
    # @4_1_P5 Waiting for a handshake response
    #
    # Covered by design, as we only get the network socket if and only if the handshake is done.

    if fix_small_message_latency
        tcpnodelay(result.io)
    end

    connection = client.connection

    io = maybewrapio(result.io)

    # For `writer` the target object is the IO stream for the WebSocket connection.
    writer = WriterProxy(io)

    # Requirement
    # @4_1_8 Handshake response is valid
    state_open(handler)

    # This function stops all the task proxies, effectively cleaning up the WSClient. This is
    # necessary when one wants to reconnect.
    cleanup = () -> begin
        stopproxy(writer)
        stopproxy(handler)
        stopproxy(connection.logic_proxy)
        stop(connection.pinger)
        if connection.reader != nothing
            stop(connection.reader)
        end
    end

    # `ClientProtocol` starts in the `STATE_OPEN` state, because it isn't responsible for making
    # connections. The target object for `logic_proxy` is the `ClientProtocol` object created here.
    framewriter = FrameWriter(writer, connection.rng)
    logic = ClientProtocol(handler, framewriter, connection.ponger, cleanup)
    connection.logic_proxy = ClientProtocolProxy(logic)

    # `Ponger` requires a logic object it can alert when a pong request hasn't been received within
    # the expected time frame. This attaches that logic object to the ponger.
    attach(connection.ponger, connection.logic_proxy)

    # `Pinger` sends ping requests at regular intervals.
    attach(connection.pinger, connection.logic_proxy)

    # The target for `reader` is the same stream we're writing to.
    connection.reader = start_reader(io, connection.logic_proxy)
    true
end

"The HTTP Upgrade failed, for whatever reason."
function connection_result_(client::WSClient,
                            result::BadHandshake,
                            handler::WebSocketHandler,
                            fix_small_message_latency::Bool)
    # Requirement
    # @4_1_EstablishConnection_4   Could not open the connection
    # @4_1_EstablishConnection_5-2 TLS Connection fails

    # Calling `state_closed` here, because that's where we're expected to attempt to reconnect.
    state_closed(handler)
    false
end

"""
Connect the client to a WebSocket server at `uri`, and use `handler` for the callbacks.

# Arguments
- `fix_small_message_latency::Bool = false`: Set the TCP_NODELAY flag to improve small message latency.

# Fix small message latency
The TCP protocol can buffer small messages (1448 bytes and smaller). The reason is that this reduces
the overhead when sending large amounts of small packets. However, it also means that latency can be
much higher for small messages. This buffering can be disabled by setting a flag TCP_NODELAY.
By default, the WebSocket client will now set the TCP_NODELAY flag.

If your application will send and receive primarily small messages (1448 bytes or smaller), and it
is sensitive to latency, then leave `fix_small_message_latency` set to true (now the default).
This sets the TCP_NODELAY flag. If you are not concerned about latency, but concerned about
throughput for many small messages, then you can set `fix_small_message_latency = false`. Then you
may get higher throughput, at the expense of higher latency for small messages.
"""
function wsconnect(client::WSClient, uri::String, handler::WebSocketHandler;
                   fix_small_message_latency=true)
    handler_proxy = WebSocketsHandlerProxy(handler)

    # Requirement
    # @4_1_P1 Initial connection state

    client.connection = WebSocketConnectionImpl()

    # The first state is always Connecting.
    state_connecting(handler_proxy, client.connection)

    # Requirement
    # @4_1_EstablishConnection_3-2 Not using a proxy

    # This makes a HTTP request to the URI and attempts to upgrade the connection to the WebSocket
    # protocol.
    handshake_result = performhandshake(client.handshake, uri)
    connection_result_(client, handshake_result, handler_proxy, fix_small_message_latency)
end

# Requirement
# @7_3-2 Clients should not close the WebSocket connection arbitrarily

"Close the WebSocket connection."
function stop(connection::WebSocketConnection)
    if connection.logic_proxy != nothing
        handle(connection.logic_proxy, CloseRequest())
    end
end

# Requirement
# @6_1-5 Opcode in the first frame
#
# Covered by design, by `send_text` and `send_binary`.

"Send a single text frame."
send_text(connection::WebSocketConnection, s::String) = handle(connection.logic_proxy, SendTextFrame(s, true, OPCODE_TEXT))

"Send a single binary frame."
send_binary(connection::WebSocketConnection, data::AbstractVector{UInt8}) =
    handle(connection.logic_proxy, SendBinaryFrame(data, true, OPCODE_BINARY))

sendmultiframetext(connection::WebSocketConnection) = TextFrameSender(connection.logic_proxy)
sendmultiframebinary(connection::WebSocketConnection) = BinaryFrameSender(connection.logic_proxy)


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: kafkaconnect
using AWS.Compat
using AWS.UUIDs

"""
    create_connector(capacity, connector_configuration, connector_name, kafka_cluster, kafka_cluster_client_authentication, kafka_cluster_encryption_in_transit, kafka_connect_version, plugins, service_execution_role_arn)
    create_connector(capacity, connector_configuration, connector_name, kafka_cluster, kafka_cluster_client_authentication, kafka_cluster_encryption_in_transit, kafka_connect_version, plugins, service_execution_role_arn, params::Dict{String,<:Any})

Creates a connector using the specified properties.

# Arguments
- `capacity`: Information about the capacity allocated to the connector. Exactly one of the
  two properties must be specified.
- `connector_configuration`: A map of keys to values that represent the configuration for
  the connector.
- `connector_name`: The name of the connector.
- `kafka_cluster`: Specifies which Apache Kafka cluster to connect to.
- `kafka_cluster_client_authentication`: Details of the client authentication used by the
  Apache Kafka cluster.
- `kafka_cluster_encryption_in_transit`: Details of encryption in transit to the Apache
  Kafka cluster.
- `kafka_connect_version`: The version of Kafka Connect. It has to be compatible with both
  the Apache Kafka cluster's version and the plugins.
- `plugins`: Specifies which plugins to use for the connector.
- `service_execution_role_arn`: The Amazon Resource Name (ARN) of the IAM role used by the
  connector to access the Amazon Web Services resources that it needs. The types of resources
  depends on the logic of the connector. For example, a connector that has Amazon S3 as a
  destination must have permissions that allow it to write to the S3 destination bucket.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"connectorDescription"`: A summary description of the connector.
- `"logDelivery"`: Details about log delivery.
- `"workerConfiguration"`: Specifies which worker configuration to use with the connector.
"""
function create_connector(
    capacity,
    connectorConfiguration,
    connectorName,
    kafkaCluster,
    kafkaClusterClientAuthentication,
    kafkaClusterEncryptionInTransit,
    kafkaConnectVersion,
    plugins,
    serviceExecutionRoleArn;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "POST",
        "/v1/connectors",
        Dict{String,Any}(
            "capacity" => capacity,
            "connectorConfiguration" => connectorConfiguration,
            "connectorName" => connectorName,
            "kafkaCluster" => kafkaCluster,
            "kafkaClusterClientAuthentication" => kafkaClusterClientAuthentication,
            "kafkaClusterEncryptionInTransit" => kafkaClusterEncryptionInTransit,
            "kafkaConnectVersion" => kafkaConnectVersion,
            "plugins" => plugins,
            "serviceExecutionRoleArn" => serviceExecutionRoleArn,
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_connector(
    capacity,
    connectorConfiguration,
    connectorName,
    kafkaCluster,
    kafkaClusterClientAuthentication,
    kafkaClusterEncryptionInTransit,
    kafkaConnectVersion,
    plugins,
    serviceExecutionRoleArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "POST",
        "/v1/connectors",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "capacity" => capacity,
                    "connectorConfiguration" => connectorConfiguration,
                    "connectorName" => connectorName,
                    "kafkaCluster" => kafkaCluster,
                    "kafkaClusterClientAuthentication" => kafkaClusterClientAuthentication,
                    "kafkaClusterEncryptionInTransit" => kafkaClusterEncryptionInTransit,
                    "kafkaConnectVersion" => kafkaConnectVersion,
                    "plugins" => plugins,
                    "serviceExecutionRoleArn" => serviceExecutionRoleArn,
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_custom_plugin(content_type, location, name)
    create_custom_plugin(content_type, location, name, params::Dict{String,<:Any})

Creates a custom plugin using the specified properties.

# Arguments
- `content_type`: The type of the plugin file.
- `location`: Information about the location of a custom plugin.
- `name`: The name of the custom plugin.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: A summary description of the custom plugin.
"""
function create_custom_plugin(
    contentType, location, name; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "POST",
        "/v1/custom-plugins",
        Dict{String,Any}(
            "contentType" => contentType, "location" => location, "name" => name
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_custom_plugin(
    contentType,
    location,
    name,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "POST",
        "/v1/custom-plugins",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "contentType" => contentType, "location" => location, "name" => name
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    create_worker_configuration(name, properties_file_content)
    create_worker_configuration(name, properties_file_content, params::Dict{String,<:Any})

Creates a worker configuration using the specified properties.

# Arguments
- `name`: The name of the worker configuration.
- `properties_file_content`: Base64 encoded contents of connect-distributed.properties file.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"description"`: A summary description of the worker configuration.
"""
function create_worker_configuration(
    name, propertiesFileContent; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "POST",
        "/v1/worker-configurations",
        Dict{String,Any}("name" => name, "propertiesFileContent" => propertiesFileContent);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_worker_configuration(
    name,
    propertiesFileContent,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "POST",
        "/v1/worker-configurations",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "name" => name, "propertiesFileContent" => propertiesFileContent
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_connector(connector_arn)
    delete_connector(connector_arn, params::Dict{String,<:Any})

Deletes the specified connector.

# Arguments
- `connector_arn`: The Amazon Resource Name (ARN) of the connector that you want to delete.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"currentVersion"`: The current version of the connector that you want to delete.
"""
function delete_connector(connectorArn; aws_config::AbstractAWSConfig=global_aws_config())
    return kafkaconnect(
        "DELETE",
        "/v1/connectors/$(connectorArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_connector(
    connectorArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "DELETE",
        "/v1/connectors/$(connectorArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_connector(connector_arn)
    describe_connector(connector_arn, params::Dict{String,<:Any})

Returns summary information about the connector.

# Arguments
- `connector_arn`: The Amazon Resource Name (ARN) of the connector that you want to
  describe.

"""
function describe_connector(connectorArn; aws_config::AbstractAWSConfig=global_aws_config())
    return kafkaconnect(
        "GET",
        "/v1/connectors/$(connectorArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_connector(
    connectorArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "GET",
        "/v1/connectors/$(connectorArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_custom_plugin(custom_plugin_arn)
    describe_custom_plugin(custom_plugin_arn, params::Dict{String,<:Any})

A summary description of the custom plugin.

# Arguments
- `custom_plugin_arn`: Returns information about a custom plugin.

"""
function describe_custom_plugin(
    customPluginArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "GET",
        "/v1/custom-plugins/$(customPluginArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_custom_plugin(
    customPluginArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "GET",
        "/v1/custom-plugins/$(customPluginArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_worker_configuration(worker_configuration_arn)
    describe_worker_configuration(worker_configuration_arn, params::Dict{String,<:Any})

Returns information about a worker configuration.

# Arguments
- `worker_configuration_arn`: The Amazon Resource Name (ARN) of the worker configuration
  that you want to get information about.

"""
function describe_worker_configuration(
    workerConfigurationArn; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "GET",
        "/v1/worker-configurations/$(workerConfigurationArn)";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function describe_worker_configuration(
    workerConfigurationArn,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "GET",
        "/v1/worker-configurations/$(workerConfigurationArn)",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_connectors()
    list_connectors(params::Dict{String,<:Any})

Returns a list of all the connectors in this account and Region. The list is limited to
connectors whose name starts with the specified prefix. The response also includes a
description of each of the listed connectors.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"connectorNamePrefix"`: The name prefix that you want to use to search for and list
  connectors.
- `"maxResults"`: The maximum number of connectors to list in one response.
- `"nextToken"`: If the response of a ListConnectors operation is truncated, it will
  include a NextToken. Send this NextToken in a subsequent request to continue listing from
  where the previous operation left off.
"""
function list_connectors(; aws_config::AbstractAWSConfig=global_aws_config())
    return kafkaconnect(
        "GET", "/v1/connectors"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_connectors(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "GET",
        "/v1/connectors",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_custom_plugins()
    list_custom_plugins(params::Dict{String,<:Any})

Returns a list of all of the custom plugins in this account and Region.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of custom plugins to list in one response.
- `"nextToken"`: If the response of a ListCustomPlugins operation is truncated, it will
  include a NextToken. Send this NextToken in a subsequent request to continue listing from
  where the previous operation left off.
"""
function list_custom_plugins(; aws_config::AbstractAWSConfig=global_aws_config())
    return kafkaconnect(
        "GET", "/v1/custom-plugins"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_custom_plugins(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "GET",
        "/v1/custom-plugins",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_worker_configurations()
    list_worker_configurations(params::Dict{String,<:Any})

Returns a list of all of the worker configurations in this account and Region.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: The maximum number of worker configurations to list in one response.
- `"nextToken"`: If the response of a ListWorkerConfigurations operation is truncated, it
  will include a NextToken. Send this NextToken in a subsequent request to continue listing
  from where the previous operation left off.
"""
function list_worker_configurations(; aws_config::AbstractAWSConfig=global_aws_config())
    return kafkaconnect(
        "GET",
        "/v1/worker-configurations";
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_worker_configurations(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return kafkaconnect(
        "GET",
        "/v1/worker-configurations",
        params;
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    update_connector(capacity, connector_arn, current_version)
    update_connector(capacity, connector_arn, current_version, params::Dict{String,<:Any})

Updates the specified connector.

# Arguments
- `capacity`: The target capacity.
- `connector_arn`: The Amazon Resource Name (ARN) of the connector that you want to update.
- `current_version`: The current version of the connector that you want to update.

"""
function update_connector(
    capacity,
    connectorArn,
    currentVersion;
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "PUT",
        "/v1/connectors/$(connectorArn)",
        Dict{String,Any}("capacity" => capacity, "currentVersion" => currentVersion);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function update_connector(
    capacity,
    connectorArn,
    currentVersion,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return kafkaconnect(
        "PUT",
        "/v1/connectors/$(connectorArn)",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "capacity" => capacity, "currentVersion" => currentVersion
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
# This elixir and indicator is only for testing purposes and does not have any practical use

using OrdinaryDiffEq
using Trixi

# Define new structs inside a module to allow re-evaluating the file.
# This module name needs to be unique among all examples, otherwise Julia will throw warnings
# if multiple test cases using the same module name are run in the same session.
module TrixiExtensionEulerAMR

using Trixi

struct IndicatorRefineCoarsen{Cache<:NamedTuple} <: Trixi.AbstractIndicator
  cache::Cache
end

function IndicatorRefineCoarsen(semi)
  basis = semi.solver.basis
  alpha = Vector{real(basis)}()
  cache = (; semi.mesh, alpha)

  return IndicatorRefineCoarsen{typeof(cache)}(cache)
end

function (indicator::IndicatorRefineCoarsen)(u::AbstractArray{<:Any,4},
                                             equations, dg, cache;
                                             t, kwargs...)
  alpha = indicator.cache.alpha
  resize!(alpha, nelements(dg, cache))

  if t >= 0.7 && t < 1.0
    # Refine to max level
    alpha .= 1.0
  elseif t >= 1.0
    # Coarsen to base level
    alpha .= -1.0
  else
    alpha .= 0.0
  end

  return alpha
end

end # module TrixiExtensionEulerAMR

import .TrixiExtensionEulerAMR

###############################################################################
# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4)

initial_condition = initial_condition_convergence_test
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

coordinates_min = (0, 0)
coordinates_max = (2, 2)
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=3,
                n_cells_max=10_000)

semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
                                    source_terms=source_terms_convergence_test)

###############################################################################
# ODE solvers, callbacks etc.

tspan = (0.0, 2.0)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

amr_controller = ControllerThreeLevel(semi, TrixiExtensionEulerAMR.IndicatorRefineCoarsen(semi),
                                      base_level=3, max_level=6,
                                      med_threshold=0.1, max_threshold=0.6)

amr_callback = AMRCallback(semi, amr_controller,
                           interval=5,
                           adapt_initial_condition=true,
                           adapt_initial_condition_only_refine=true)

stepsize_callback = StepsizeCallback(cfl=1.0)

callbacks = CallbackSet(summary_callback,
                        analysis_callback, alive_callback,
                        save_solution, amr_callback,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary