text:
using KernelAbstractions: NoneEvent
using OffsetArrays: OffsetArray

using Oceananigans.Operators: Δzᵃᵃᶜ
using Oceananigans.BoundaryConditions: left_gradient, right_gradient, linearly_extrapolate, FBC, VBC, GBC
using Oceananigans.BoundaryConditions: fill_bottom_halo!, fill_top_halo!, apply_z_bottom_bc!, apply_z_top_bc!
using Oceananigans.Grids: Flat, Bounded
using Oceananigans.Architectures: device_event

import Oceananigans.Utils: launch!
import Oceananigans.Grids: validate_size, validate_halo
import Oceananigans.BoundaryConditions: fill_halo_regions!
import Oceananigans.TurbulenceClosures: time_discretization
import Oceananigans.TurbulenceClosures: calculate_diffusivities!

#####
##### Implements a "single column model mode" for HydrostaticFreeSurfaceModel
#####

const SingleColumnGrid = AbstractGrid{<:AbstractFloat, <:Flat, <:Flat, <:Bounded}

#####
##### Model constructor utils
#####

PressureField(arch, ::SingleColumnGrid) = (pHY′ = nothing,)
FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, velocities, arch, ::SingleColumnGrid) = nothing
FreeSurface(free_surface::ImplicitFreeSurface{Nothing}, velocities, arch, ::SingleColumnGrid) = nothing

validate_momentum_advection(momentum_advection, ::SingleColumnGrid) = nothing
validate_tracer_advection(tracer_advection::AbstractAdvectionScheme, ::SingleColumnGrid) = nothing, NamedTuple()
validate_tracer_advection(tracer_advection::Nothing, ::SingleColumnGrid) = nothing, NamedTuple()

#####
##### Time-step optimizations
#####

calculate_free_surface_tendency!(arch, ::SingleColumnGrid, args...) = NoneEvent()

# Fast state update and halo filling

function update_state!(model::HydrostaticFreeSurfaceModel, grid::SingleColumnGrid)

    fill_halo_regions!(prognostic_fields(model), model.architecture, model.clock, fields(model))

    compute_auxiliary_fields!(model.auxiliary_fields)

    # Calculate diffusivities
    calculate_diffusivities!(model.diffusivity_fields, model.closure, model)

    fill_halo_regions!(model.diffusivity_fields,
                       model.architecture,
                       model.clock,
                       fields(model))

    return nothing
end

import Oceananigans.TurbulenceClosures: ∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ, ∂ⱼ_τ₃ⱼ, ∇_dot_qᶜ

using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure

const ClosureArray = AbstractArray{<:AbstractTurbulenceClosure}

@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid::SingleColumnGrid, closure_array::ClosureArray, args...)
    @inbounds closure = closure_array[i, j]
    return ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, args...)
end

@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid::SingleColumnGrid, closure_array::ClosureArray, args...)
    @inbounds closure = closure_array[i, j]
    return ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure, args...)
end

@inline function ∇_dot_qᶜ(i, j, k, grid::SingleColumnGrid, closure_array::ClosureArray, c, tracer_index, args...)
    @inbounds closure = closure_array[i, j]
    return ∇_dot_qᶜ(i, j, k, grid, closure, c, tracer_index, args...)
end

struct ColumnEnsembleSize{C<:Tuple{Int, Int}}
    ensemble :: C
    Nz :: Int
    Hz :: Int
end

ColumnEnsembleSize(; Nz, ensemble=(0, 0), Hz=1) = ColumnEnsembleSize(ensemble, Nz, Hz)

validate_size(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(e.ensemble[1], e.ensemble[2], e.Nz)
validate_halo(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(0, 0, e.Hz)

@inline time_discretization(::AbstractArray{<:AbstractTurbulenceClosure{TD}}) where TD = TD()

#####
##### TKEBasedVerticalDiffusivity helpers
#####

using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVD, Kuᶜᶜᶜ, Kcᶜᶜᶜ, Keᶜᶜᶜ, _top_tke_flux, CATKEVDArray

import Oceananigans.TurbulenceClosures: calculate_diffusivities!
import Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: top_tke_flux

function calculate_diffusivities!(diffusivities, closure::CATKEVDArray, model)

    arch = model.architecture
    grid = model.grid
    velocities = model.velocities
    tracers = model.tracers
    buoyancy = model.buoyancy
    clock = model.clock
    top_tracer_bcs = NamedTuple(c => tracers[c].boundary_conditions.top for c in propertynames(tracers))

    event = launch!(arch, grid, :xyz,
                    calculate_CATKEArray_diffusivities!,
                    diffusivities, grid, closure, velocities, tracers, buoyancy, clock, top_tracer_bcs,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    return nothing
end

@kernel function calculate_CATKEArray_diffusivities!(diffusivities, grid::SingleColumnGrid, closure_array::CATKEVDArray, args...)
    i, j, k, = @index(Global, NTuple)
    @inbounds begin
        closure = closure_array[i, j]
        diffusivities.Kᵘ[i, j, k] = Kuᶜᶜᶜ(i, j, k, grid, closure, args...)
        diffusivities.Kᶜ[i, j, k] = Kcᶜᶜᶜ(i, j, k, grid, closure, args...)
        diffusivities.Kᵉ[i, j, k] = Keᶜᶜᶜ(i, j, k, grid, closure, args...)
    end
end

@inline tracer_tendency_kernel_function(model::HydrostaticFreeSurfaceModel, closure::CATKEVDArray, ::Val{:e}) =
    hydrostatic_turbulent_kinetic_energy_tendency

""" Compute the flux of TKE through the surface / top boundary. """
@inline function top_tke_flux(i, j, grid::SingleColumnGrid, clock, fields, parameters, closure_array::CATKEVDArray, buoyancy)
    top_tracer_bcs = parameters.top_tracer_boundary_conditions
    top_velocity_bcs = parameters.top_velocity_boundary_conditions
    @inbounds closure = closure_array[i, j]

    return _top_tke_flux(i, j, grid, closure.surface_TKE_flux, closure,
                         buoyancy, fields, top_tracer_bcs, top_velocity_bcs, clock)
end

@inline function hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid::SingleColumnGrid,
                                                               val_tracer_index::Val{tracer_index},
                                                               advection,
                                                               closure_array::CATKEVDArray, args...) where tracer_index

    @inbounds closure = closure_array[i, j]
    return hydrostatic_turbulent_kinetic_energy_tendency(i, j, k, grid, val_tracer_index, advection, closure, args...)
end

#####
##### Arrays of Coriolises
#####

import Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U

const FPlaneArray = AbstractArray{<:FPlane, 2}

@inline x_f_cross_U(i, j, k, grid::SingleColumnGrid, coriolis::FPlaneArray, U) = @inbounds - coriolis[i, j].f * ℑxyᶠᶜᵃ(i, j, k, grid, U[2])
@inline y_f_cross_U(i, j, k, grid::SingleColumnGrid, coriolis::FPlaneArray, U) = @inbounds   coriolis[i, j].f * ℑxyᶜᶠᵃ(i, j, k, grid, U[1])
@inline z_f_cross_U(i, j, k, grid::SingleColumnGrid, coriolis::FPlaneArray, U) = zero(eltype(grid))


text:
using Test
using Statistics
using BenchmarkTools

using Cyclotomics

function zumbroich_perf(ran)
    @testset "zumbroich perf" begin
        @test all(
            Cyclotomics.zumbroich_plain.(ran) .==
            sort!.(collect.(Cyclotomics.zumbroich_basis.(ran))) .==
            Cyclotomics.zumbroich_direct.(ran),
        )

        plain = [(@timed sum(first, Cyclotomics.zumbroich_plain(i)))[2] for i in ran]
        @info "By plain loop: mean and variance over $ran" μ =
            round(mean(plain), sigdigits = 4) σ =
            round(Statistics.var(plain), sigdigits = 4)
        @btime Cyclotomics.zumbroich_plain($(last(ran)))

        compl = [(@timed sum(first, Cyclotomics.zumbroich_basis(i)))[2] for i in ran]
        @info "Constructing the complement: mean and variance over $ran" μ =
            round(mean(compl), sigdigits = 4) σ =
            round(Statistics.var(compl), sigdigits = 4)
        @btime Cyclotomics.zumbroich_basis($(last(ran)))

        direc = [(@timed length(Cyclotomics.zumbroich_direct(i)))[2] for i in ran]
        @info "GAP style: mean and variance over $ran" μ = round(mean(direc), sigdigits = 4) σ =
            round(Statistics.var(direc), sigdigits = 4)
        @btime Cyclotomics.zumbroich_direct($(last(ran)))
    end
end

zumbroich_perf(1:1000)
zumbroich_perf(25000:30001)

ext: jl


text:
@testset "Ensembles" begin
  d = PointSet(rand(2,100))
  r = (z=[1:100 for i in 1:10],)
  s = Ensemble(d, r)
  @test domain(s) == d
  @test s[:z] == r[:z]
  for i in 1:10
    @test s[i] == georef((z=1:100,), d)
  end

  @test sprint(show, s) == "2D Ensemble"
  @test sprint(show, MIME"text/plain"(), s) == "2D Ensemble\n  domain: 100 PointSet{2,Float64}\n  variables: z\n  N° reals:  10"

  d = CartesianGrid(10,10)
  r = (z=[1:100 for i in 1:10],)
  s = Ensemble(d, r)
  @test domain(s) == d
  @test s[:z] == r[:z]
  for i in 1:10
    @test s[i] == georef((z=1:100,), d)
  end

  @test sprint(show, s) == "2D Ensemble"
  @test sprint(show, MIME"text/plain"(), s) == "2D Ensemble\n  domain: 10×10 CartesianGrid{2,Float64}\n  variables: z\n  N° reals:  10"

  if visualtests
    @test_reference "data/ensemble.png" plot(s,size=(800,300))
  end
end


text:
#clearconsole()
push!(LOAD_PATH,"../src")
using Pkg
using Dates
using DataFrames
using Documenter, DocumenterTools
using Distributions, Statistics, StatsBase
using MCHammer
using Gadfly, Compose, Cairo

makedocs(
sitename="MCHammer.jl",
modules =[MCHammer],
pages = Any[
"Home" => "index.md",
"User Manual" => Any[
    "Simulation Functions" => "manual/1_functions.md",
    "Charting & Analyzing" => "manual/2_charts.md",
    "Time-Series Functions" => "manual/3_time_series.md",
],
"Tutorials" => Any[
    "My First Model" => "tutorials/1_first_model.md",
    "Correlating Inputs" => "tutorials/2_Correlated_Model.md",
    "Simulated CashFlow Model" => "tutorials/3_NPV_testmodel.md",
],
"Import / Export Data" => "manual/4_moving_results.md"
],

format = Documenter.HTML(
    # Use clean URLs, unless built as a "local" build
    prettyurls = get(ENV, "CI", nothing) == "true",
    canonical = "https://etorkia.github.io/MCHammer.jl/",
    assets = ["assets/favicon.ico"],
    analytics = "UA-3913053-5",
),

authors = "Eric Torkia, Technology Partnerz and contributors",
doctest = true,
source  = "src",
build   = "build",
clean   = true,
repo = "github.com/etorkia/MCHammer.jl")

deploydocs(
    target="build",
    repo = "github.com/etorkia/MCHammer.jl.git",
    branch = "gh-pages",
    devbranch="master"
    versions = ["stable" => "v^", "v#.#"]
)

#deploydocs(repo = "github.com/etorkia/MCHammer.jl.git")


text:
using WFDB
using Test

@testset "WFDB.jl" begin
    # Write your tests here.
end


text:
import Base: Float64, Float32, Float16, Int

primitive type FastFloat16 <: AbstractFloat 32 end

# CONVERSIONS
# helper functions to create correct shift & set/shavemasks
shift(n::UInt32) = 13 + n                   # for subnormals n > 0 to decrease precision
mask(n::UInt32) = 0x003f_ffff >> (0xa-n)    # same: mask more significant bits

# conversion FastFloat16 -> Float32 is just reinterpretation
Base.Float32(x::FastFloat16) = reinterpret(Float32,x)

# conversion Float32 -> via rounding & over/underflow cap
function FastFloat16(x::Float32)
    ui = reinterpret(UInt32,x)
    sign = ui & 0x8000_0000     # separate the sign bit
    ui = ui & 0x7fff_ffff       # from the exp&significant bits

    # check for Float16-subnormals: subnormal sn = ui < 0x3880_0000
    scale = 0x71 - (ui >> 23)           # 0 for sn, 1 for sn/2, 2 for sn/4 etc.
    # for x>=subnormal mask&shift operators are constant, for subnormals the
    # precision is decreased by 1 bit for every power 2
    maskbits,shiftbits = scale < 0xb ?  # = subnormal? yes: get variable mask&shift, no: constant
        (mask(scale), shift(scale)) : (0x0000_0fff, 13)

    # round to nearest (retain 10 sigbits for >subnormal, less for subnormals)
    ui += maskbits + ((ui >> shiftbits) & 0x0000_0001)  # take carry over into account
    ui &= (~maskbits) << 1          # set the non-Float16 sigbits to 0

    # check for overflow/underflow
    # underflow = ui smaller than minpos/2 = Float32(nextfloat(Float16(0)))/2
    ui = ui < 0x3300_0000 ? 0x0000_0000 : ui   # then replace with 0
    ui = ui > 0x477f_e000 ? 0x7f80_0000 : ui   # overflow? replace with Inf32

    # recombine sign&exp&significant
    return reinterpret(FastFloat16,sign | ui)
end

# conversion between FastFloatt16 and various floats
Float64(x::FastFloat16) = Float64(Float32(x))
Float16(x::FastFloat16) = Float16(Float32(x))
FastFloat16(x::Float64) = FastFloat16(Float32(x))
FastFloat16(x::Float16) = FastFloat16(Float32(x))
Int(x::FastFloat16) = Int(Float32(x))
FastFloat16(x::Int) = FastFloat16(Float32(x))

Base.iszero(x::FastFloat16) = iszero(Float32(x))
Base.isnan(x::FastFloat16) = isnan(Float32(x))
Base.signbit(x::FastFloat16) = signbit(Float32(x))

Base.zero(::Type{FastFloat16}) = reinterpret(FastFloat16,0x0000_0000)
Base.one(::Type{FastFloat16}) = reinterpret(FastFloat16,0x3f80_0000)

Base.floatmin(::Type{FastFloat16}) = reinterpret(FastFloat16,0x3880_0000)   # = Float32(floatmin(Float16))
Base.floatmax(::Type{FastFloat16}) = reinterpret(FastFloat16,0x477f_e000)   # = Float32(floatmax(Float16))

# In the absence of -Inf,Inf define typemin,typemax as floatmin,floatmax
Base.typemin(::Type{FastFloat16}) = reinterpret(FastFloat16,0xff800000)     # = Float32(typemin(Float16))
Base.typemax(::Type{FastFloat16}) = reinterpret(FastFloat16,0x7f800000)     # = Float32(typemax(Float16))

Base.:(-)(x::FastFloat16) = FastFloat16(-Float32(x))
Base.inv(x::FastFloat16) = FastFloat16(inv(Float32(x)))

# Arithmetics
Base.:(*)(x::FastFloat16,y::FastFloat16) = FastFloat16(Float32(x)*Float32(y))
Base.:(/)(x::FastFloat16,y::FastFloat16) = FastFloat16(Float32(x)/Float32(y))
Base.:(+)(x::FastFloat16,y::FastFloat16) = FastFloat16(Float32(x)+Float32(y))
Base.:(-)(x::FastFloat16,y::FastFloat16) = FastFloat16(Float32(x)-Float32(y))
Base.sqrt(x::FastFloat16) = FastFloat16(sqrt(Float32(x)))

Base.round(x::FastFloat16, r::RoundingMode{:Up}) = Int(ceil(Float32(x)))
Base.round(x::FastFloat16, r::RoundingMode{:Down}) = Int(floor(Float32(x)))
Base.round(x::FastFloat16, r::RoundingMode{:Nearest}) = Int(round(Float32(x)))

for t in (Int8, Int16, Int32, Int64, Int128, UInt8, UInt16, UInt32, UInt64, UInt128)
	@eval Base.promote_rule(::Type{FastFloat16}, ::Type{$t}) = FastFloat16
end

Base.nextfloat(x::FastFloat16) = FastFloat16(nextfloat(Float16(x)))
Base.nextfloat(x::FastFloat16,n::Int) = FastFloat16(nextfloat(Float16(x),n))
Base.prevfloat(x::FastFloat16) = FastFloat16(prevfloat(Float16(x)))
Base.prevfloat(x::FastFloat16,n::Int) = FastFloat16(prevfloat(Float16(x),n))

Base.log2(x::FastFloat16) = FastFloat16(log2(Float32(x)))
Base.log(x::FastFloat16) = FastFloat16(log(Float32(x)))

Base.:(==)(x::FastFloat16,y::FastFloat16) = Float32(x) == Float32(y)
Base.:(>)(x::FastFloat16,y::FastFloat16) = Float32(x) > Float32(y)
Base.:(<)(x::FastFloat16,y::FastFloat16) = Float32(x) < Float32(y)
Base.:(<=)(x::FastFloat16,y::FastFloat16) = Float32(x) <= Float32(y)

# Showing
function Base.show(io::IO, x::FastFloat16)
    io2 = IOBuffer()
    print(io2,Float16(x))
    f = String(take!(io2))
    print(io,"FastFloat16("*f*")")
end

Base.bitstring(x::FastFloat16) = bitstring(Float32(x))

function Base.bitstring(x::FastFloat16,mode::Symbol)
    if mode == :split	# split into sign, integer, fraction
        s = bitstring(x)
		return "$(s[1]) $(s[2:9]) $(s[10:end])"
    else
        return bitstring(x)
    end
end


text:
using Nanoflann
using Base.Test

function bench(tree, query_points)
    knn(tree, query_points, 10)
end

N = 100000
L = 10
data = L*rand(Float64, 3, N)
query_points = L*rand(Float64, 3, N)
tree1 = Nanoflann.KDTree(data)
tree2 = NearestNeighbors.KDTree(data)

@time Nanoflann.KDTree(data)
@time NearestNeighbors.KDTree(data)

bench(tree1, query_points[:,1:2])
bench(tree2, query_points[:,1:2])
@time inds1,dist1 = bench(tree1, query_points)
@time inds2,dist2 = bench(tree2, query_points)

for i = 1:N
    @assert vec(inds1[sortperm(dist1[:,i]),i]) == inds2[i][sortperm(dist2[i])]
end

# using Displaz
# inds,_ = knn(tree1, 0.5*L*ones(3), 100)
# Displaz.clf()
# Displaz.plot(data')
# Displaz.hold(true)
# Displaz.plot([0.5,0.5,0.5]', color="b")
# Displaz.plot(data[:,inds]', color="r")

nothing


text:
s = """
module test; if true; end
"""
msgs = lintstr(s)
@test msgs[1].code == :E112
@test contains(msgs[1].message, "incomplete:")

s = """
2 +
"""
msgs = lintstr(s)
@test msgs[1].code == :E112
@test contains(msgs[1].message, "incomplete:")


text:
const BUFFER_TYPE = @map {Symbol, GLenum} {
    Array       = GL_ARRAY_BUFFER,
    Element     = GL_ELEMENT_ARRAY_BUFFER,
    Texture     = GL_TEXTURE_BUFFER,
    Uniform     = GL_UNIFORM_BUFFER,
    Atomic      = GL_ATOMIC_COUNTER_BUFFER,
    CopyRead    = GL_COPY_READ_BUFFER,
    CopyWrite   = GL_COPY_WRITE_BUFFER,
    Dispatch    = GL_DISPATCH_INDIRECT_BUFFER,
    Draw        = GL_DRAW_INDIRECT_BUFFER,
    PixelPack   = GL_PIXEL_PACK_BUFFER,
    PixelUnpack = GL_PIXEL_UNPACK_BUFFER,
}
"""
        genbuffer() --> GLuint

Creates a OpenGL `buffer` object.
"""
function genbuffer()
    buffer = GLuint'[]
    glGenBuffers(1, buffer)
    return buffer'
end
"""
        bindbuffer(buffer, type::Symbol) -> Nothing

Bind the `buffer` object.

The `type` parameter can be:

- `:Array`       for `GL_ARRAY_BUFFER`
- `:Element`     for `GL_ELEMENT_ARRAY_BUFFER`
- `:Texture`     for `GL_TEXTURE_BUFFER`
- `:Uniform`     for `GL_UNIFORM_BUFFER`
- `:Atomic`      for `GL_ATOMIC_COUNTER_BUFFER`
- `:CopyRead`    for `GL_COPY_READ_BUFFER`
- `:CopyWrite`   for `GL_COPY_WRITE_BUFFER`
- `:Dispatch`    for `GL_DISPATCH_INDIRECT_BUFFER`
- `:Draw`        for `GL_DRAW_INDIRECT_BUFFER`
- `:PixelPack`   for `GL_PIXEL_PACK_BUFFER`
- `:PixelUnpack` for `GL_PIXEL_UNPACK_BUFFER`
"""
function bindbuffer(type::Symbol, buffer::GLuint)
    @assert type in keys(BUFFER_TYPE) "Unknown buffer type: $type"
    glBindBuffer(BUFFER_TYPE[type], buffer)
end
"""
        bufferdata(type, data, usage) -> Nothing

Creates and initializes a `buffer` object's data store.

The `type` parameter can be:

- `:Array`       for `GL_ARRAY_BUFFER`
- `:Element`     for `GL_ELEMENT_ARRAY_BUFFER`
- `:Texture`     for `GL_TEXTURE_BUFFER`
- `:Uniform`     for `GL_UNIFORM_BUFFER`
- `:Atomic`      for `GL_ATOMIC_COUNTER_BUFFER`
- `:CopyRead`    for `GL_COPY_READ_BUFFER`
- `:CopyWrite`   for `GL_COPY_WRITE_BUFFER`
- `:Dispatch`    for `GL_DISPATCH_INDIRECT_BUFFER`
- `:Draw`        for `GL_DRAW_INDIRECT_BUFFER`
- `:PixelPack`   for `GL_PIXEL_PACK_BUFFER`
- `:PixelUnpack` for `GL_PIXEL_UNPACK_BUFFER`
"""
function bufferdata(type::Symbol, data::AbstractArray, usage::GLenum)
    target = BUFFER_TYPE[type]
    glBufferData(target, sizeof(data), data, usage)
end
"""
        arraybuffer(data; usage::Symbol) -> GLuint

Creates an `Array` buffer object and initializes with
the specified `usage` and `data`.
"""
function arraybuffer(data::AbstractArray, usage::GLenum)
    buffer = genbuffer()
    bindbuffer(:Array, buffer)
    bufferdata(:Array, data, usage)
    return buffer
end
"""
        indexbuffer(data; usage::Symbol) -> GLuint

Creates an `ElementArray` buffer object and initializes with
the specified `usage` and `data`.
"""
function indexbuffer(data::AbstractArray, usage::GLenum)
    buffer = genbuffer()
    bindbuffer(:Element, buffer)
    bufferdata(:Element, data, usage)
    return buffer
end

function vertexarray()
    id = GLuint'[0]
    glGenVertexArrays(1, id)
    glBindVertexArray(id')
    return id'
end
"""
        vertexattrib(type, stride, attributes) -> Nothing

Define an array of vertex attribute data.

- *@param* `type` - the vertex data type.
- *@param* `stride` - the size of each vertex.
- *@param* `attributes` - a map (Dict) containg the index, size and offset of each attribute.

The `attributes` parameter is a `Dict{Int, NamedTuple}`:

```julia
attributes = Dict(
    index₁ => ( size = size₁, offset = offset₁ [, normalized = false] ),
    index₂ => ( size = size₂, offset = offset₂ [, normalized = false] ),
        ⋮
    indexₙ  => ( size = sizeₙ, offset = offsetₙ  [, normalized = false] ),
)
```

# Usage

```julia
vertexattrib(Cfloat, 4, @map {
    0 := (size = 2, offset = 0),
    1 := (size = 2, offset = 2),
})
```
"""
function vertexattrib(::Type{T}, stride::Int, attributes::Dict{Int, <:NamedTuple}) where {T}
    for (index, attr) in attributes
        size = attr.size
        offset = attr.offset
        norm = hasproperty(attr, :normalized) ? attr.normalized : false
        glEnableVertexAttribArray(index)
        glVertexAttribPointer(
            index,
            size,
            genumof(T),
            gbool(norm),
            stride * sizeof(T),
            offset * sizeof(T) |> Ptr{Cvoid}
        )
    end
end

function vertexattrib(::Type{T}, stride::Int, attributes::NamedTuple...) where {T}
    indices = 0:length(attributes) - 1
    vertexattrib(T, stride, Dict(zip(indices, attributes)))
end


text:
@with_kw struct UniformSource{GeometryType,E <: AbstractExcitation} <: AbstractSource
    geometry::GeometryType
    excitation::E
    polarization::Vec
end


text:
using MPI
using StaticArrays
using ClimateMachine
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.MPIStateArrays
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using Printf
using LinearAlgebra
using Logging

using ClimateMachine.BalanceLaws

import ClimateMachine.BalanceLaws:
    vars_state,
    flux_first_order!,
    flux_second_order!,
    source!,
    wavespeed,
    boundary_state!,
    nodal_init_state_auxiliary!,
    init_state_prognostic!,
    update_auxiliary_state!,
    indefinite_stack_integral!,
    reverse_indefinite_stack_integral!,
    integral_load_auxiliary_state!,
    integral_set_auxiliary_state!,
    reverse_integral_load_auxiliary_state!,
    reverse_integral_set_auxiliary_state!

import ClimateMachine.DGMethods: init_ode_state
using ClimateMachine.Mesh.Geometry: LocalGeometry

struct IntegralTestModel{dim} <: BalanceLaw end

vars_state(::IntegralTestModel, ::DownwardIntegrals, T) = @vars(a::T, b::T)
vars_state(::IntegralTestModel, ::UpwardIntegrals, T) = @vars(a::T, b::T)
vars_state(m::IntegralTestModel, ::Auxiliary, T) = @vars(
    int::vars_state(m, UpwardIntegrals(), T),
    rev_int::vars_state(m, DownwardIntegrals(), T),
    coord::SVector{3, T},
    a::T,
    b::T,
    rev_a::T,
    rev_b::T
)

vars_state(::IntegralTestModel, ::AbstractStateType, T) = @vars()

flux_first_order!(::IntegralTestModel, _...) = nothing
flux_second_order!(::IntegralTestModel, _...) = nothing
source!(::IntegralTestModel, _...) = nothing
boundary_state!(_, ::IntegralTestModel, _...) = nothing
init_state_prognostic!(::IntegralTestModel, _...) = nothing
wavespeed(::IntegralTestModel, _...) = 1

function nodal_init_state_auxiliary!(
    ::IntegralTestModel{dim},
    aux::Vars,
    tmp::Vars,
    g::LocalGeometry,
) where {dim}
    x, y, z = aux.coord = g.coord
    if dim == 2
        aux.a = x * y
        aux.b = 2 * x * y + sin(x) * y^2 / 2 - (z - 1)^2 * y^3 / 3
        y_top = 3
        a_top = x * y_top
        b_top = 2 * x * y_top + sin(x) * y_top^2 / 2 - (z - 1)^2 * y_top^3 / 3
        aux.rev_a = a_top - aux.a
        aux.rev_b = b_top - aux.b
    else
        aux.a = x * z + y * z
        aux.b = 2 * x * z + sin(x) * y * z - (1 + (z - 1)^3) * y^2 / 3
        zz_top = 3
        a_top = x * zz_top + y * zz_top
        b_top =
            2 * x * zz_top + sin(x) * y * zz_top -
            (1 + (zz_top - 1)^3) * y^2 / 3
        aux.rev_a = a_top - aux.a
        aux.rev_b = b_top - aux.b
    end
end

function update_auxiliary_state!(
    dg::DGModel,
    m::IntegralTestModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    indefinite_stack_integral!(dg, m, Q, dg.state_auxiliary, t, elems)
    reverse_indefinite_stack_integral!(dg, m, Q, dg.state_auxiliary, t, elems)

    return true
end

@inline function integral_load_auxiliary_state!(
    m::IntegralTestModel{dim},
    integrand::Vars,
    state::Vars,
    aux::Vars,
) where {dim}
    x, y, z = aux.coord
    integrand.a = x + (dim == 3 ? y : 0)
    integrand.b = 2 * x + sin(x) * y - (z - 1)^2 * y^2
end

@inline function integral_set_auxiliary_state!(
    m::IntegralTestModel,
    aux::Vars,
    integral::Vars,
)
    aux.int.a = integral.a
    aux.int.b = integral.b
end

@inline function reverse_integral_load_auxiliary_state!(
    m::IntegralTestModel,
    integral::Vars,
    state::Vars,
    aux::Vars,
)
    integral.a = aux.int.a
    integral.b = aux.int.b
end

@inline function reverse_integral_set_auxiliary_state!(
    m::IntegralTestModel,
    aux::Vars,
    integral::Vars,
)
    aux.rev_int.a = integral.a
    aux.rev_int.b = integral.b
end

using Test
function test_run(mpicomm, dim, Ne, N, FT, ArrayType)

    brickrange = ntuple(j -> range(FT(0); length = Ne[j] + 1, stop = 3), dim)
    topl = StackedBrickTopology(
        mpicomm,
        brickrange,
        periodicity = ntuple(j -> true, dim),
    )

    grid = DiscontinuousSpectralElementGrid(
        topl,
        FloatType = FT,
        DeviceArray = ArrayType,
        polynomialorder = N,
    )
    dg = DGModel(
        IntegralTestModel{dim}(),
        grid,
        RusanovNumericalFlux(),
        CentralNumericalFluxSecondOrder(),
        CentralNumericalFluxGradient(),
    )

    Q = init_ode_state(dg, FT(0))
    dQdt = similar(Q)

    dg(dQdt, Q, nothing, 0.0)

    # Wrapping in Array ensure both GPU and CPU code use same approx
    if N[end] > 0
        # Forward integral a
        @test Array(dg.state_auxiliary.data[:, 1, :]) ≈
              Array(dg.state_auxiliary.data[:, 8, :])
        # Forward integral b
        @test Array(dg.state_auxiliary.data[:, 2, :]) ≈
              Array(dg.state_auxiliary.data[:, 9, :])
        # Reverse integral a
        @test Array(dg.state_auxiliary.data[:, 3, :]) ≈
              Array(dg.state_auxiliary.data[:, 10, :])
        # Reverse integral b
        @test Array(dg.state_auxiliary.data[:, 4, :]) ≈
              Array(dg.state_auxiliary.data[:, 11, :])
    else
        # For N = 0 we only compare the first integral which is the integral of
        # a vertical constant function; N = 0 can also integrate linears exactly
        # since we use the midpoint rule to compute the face value, but the
        # averaging procedure we use below does not work in this case
        Nq = polynomialorders(grid) .+ 1

        # Reshape the data array to be (dofs, vertical elem, horizontal elem)
        A_faces = reshape(
            Array(dg.state_auxiliary.data[:, 1, :]),
            prod(Nq),
            Ne[end],               # Vertical element is fastest on stacked meshes
            prod(Ne[1:(end - 1)]), # Horiztonal elements
        )
        A_center_exact = reshape(
            Array(dg.state_auxiliary.data[:, 8, :]),
            prod(Nq),
            Ne[end],               # Vertical element is fastest on stacked meshes
            prod(Ne[1:(end - 1)]), # Horiztonal elements
        )
        # With N = 0, the integral will return the values in the faces. Namely,
        # verical element index `eV` will be the value of the integral on the
        # face ABOVE element `eV`. Namely,
        #    A_faces[n, eV, eH]
        # will be degree of freedom `n`, in horizontal element stack `eH`, and
        # face `eV + 1/2`.
        #
        # The exact values stored in `A_center_exact` are actually at the cell
        # centers because these are computed using the `init_state_auxiliary!`
        # which evaluates using the cell centers.
        #
        # This mismatch means we need to convert from faces to cell centers for
        # comparison, and we do this using averaging to go from faces to cell
        # centers.

        # Storage for the averaging
        A_center = similar(A_faces)

        # Bottom cell value is average of 0 and top face of cell
        A_center[:, 1, :] .= A_faces[:, 1, :] / 2

        # Remaining cells are average of the two faces
        A_center[:, 2:end, :, :] .=
            (A_faces[:, 1:(end - 1), :] + A_faces[:, 2:end, :]) / 2

        # Compare the exact and computed
        @test A_center ≈ A_center_exact

        # We do the same things for the reverse integral, the only difference is
        # now the values
        #    RA_faces[n, eV, eH]
        # will be degree of freedom `n`, in horizontal element stack `eH`, and
        # face `eV - 1/2` (e.g., the face below element `eV`
        # Reshape the data array to be (dofs, vertical elm, horizontal elm)
        RA_faces = reshape(
            Array(dg.state_auxiliary.data[:, 3, :]),
            prod(Nq),
            Ne[end],               # Vertical element is fastest on stacked meshes
            prod(Ne[1:(end - 1)]), # Horiztonal elements
        )
        RA_center_exact = reshape(
            Array(dg.state_auxiliary.data[:, 10, :]),
            prod(Nq),
            Ne[end],               # Vertical element is fastest on stacked meshes
            prod(Ne[1:(end - 1)]), # Horiztonal elements
        )

        # Storage for the averaging
        RA_center = similar(RA_faces)

        # Top cell value is average of 0 and top face of cell
        RA_center[:, end, :] .= RA_faces[:, end, :] / 2

        # Remaining cells are average of the two faces
        RA_center[:, 1:(end - 1), :, :] .=
            (RA_faces[:, 1:(end - 1), :] + RA_faces[:, 2:end, :]) / 2

        # Compare the exact and computed
        @test RA_center ≈ RA_center_exact
    end
end

let
    ClimateMachine.init()
    ArrayType = ClimateMachine.array_type()

    mpicomm = MPI.COMM_WORLD

    numelem = (5, 6, 7)
    lvls = 1

    for polynomialorder in ((4, 4), (4, 3), (4, 0))
        for FT in (Float64,)
            for dim in 2:3
                err = zeros(FT, lvls)
                for l in 1:lvls
                    @info (ArrayType, FT, dim, polynomialorder)
                    test_run(
                        mpicomm,
                        dim,
                        ntuple(j -> 2^(l - 1) * numelem[j], dim),
                        polynomialorder,
                        FT,
                        ArrayType,
                    )
                end
            end
        end
    end
end

nothing

max_stars_repo_name: jm-c/CLIMA

max_issues_repo_name: jm-c/CLIMA

max_forks_repo_name: jm-c/CLIMA


text:
## docs/make.jl --- to make the docs

### got this line from DynamicalSystems.jl
CI = get(ENV, "CI", nothing) === true || !isnothing(get(ENV, "GITHUB_TOKEN", nothing))

using Documenter
using ProtStructs

makedocs(
    sitename="ProtStructs",
    authors="Bon Leif Amalla and contributors",
    doctest=false,
    format = Documenter.HTML(prettyurls=CI),
    modules = [ProtStructs,],
    pages= [
        "Introduction" => "index.md",
        "Brief walkthrough" => "start.md",
        "Models" => [
            "Structure Frames" => "structs/sfs.md",
            # "HBond Dictionary" => "structs/hbdict.md"
        ],
        # "Hydrogen bonding model" => "hbond.md",
        # "Secondary structure" => "secstruct.md",
        # "Implementation notes" => "impnotes.md",
        "API Documentation" => "API.md"
   ]
)

if CI
    deploydocs(
        repo = "https://github.com/bldamalla/ProtStructs.jl.git",
        target = "build",
        push_preview = false
    )
end


text:
struct BitStream
	data::Array{Int8, 1}
	indices::Array{Int, 1}
	BitStream() = new(Array{Int8, 1}(undef, 16), [1, 1])
	BitStream(prealloc::Int) = new(Array{Int8, 1}(undef, prealloc), [1, 1])
end

function push!(strm::BitStream, data::Bool)
	strm.data[strm.indices[1]] |= (data ? 1 : 0) << (strm.indices[2] - 1)

	if strm.indices[2] == 8
		strm.indices[2] = 1
		strm.indices[1] += 1
		push!(strm.data, 0)
	end
end

function get(strm::BitStream, index::Int)::Bool
	return (strm.data[cld(index, 8)] & (1 << ((index - 1) % 8))) != 0
end

function set(strm::BitStream, index::Int, val::Bool)
	strm.data[cld(index, 8)] &= ~(1 << ((index - 1) % 8))
	strm.data[cld(index, 8)] |= (val ? 1 : 0) << ((index - 1) % 8)
end

function length(strm::BitStream)::Int
	return (strm.indices[1] - 1) * 8 + strm.indices[2]
end

function print(strm::BitStream)
	for i in 1 : length(strm)
		print(get(strm, i) ? '1' : 0)
	end
end

function println(strm::BitStream)
	print(strm)
	print('\n')
end


text:
"Meta data to attach to ω::Ω"
const Tags{K, V} = NamedTuple{K, V}

"Does tag type contain `t`, forall t in `tags`?"
hastags(::Type{Tags{K, V}}, tags::Symbol...) where {K, V} = all([t in K for t in tags])
combinetag(::Type{Val{:replmap}}, a, b) = merge(a, b)

"Sample space tagged with meta-data.  Enables `iid`, `replace`, `trackerror`"
struct TaggedΩ{I, TAGS <: Tags, ΩT} <: Ω{I}
  taggedω::ΩT
  tags::TAGS
  TaggedΩ(ω::ΩT, tags::TAGS) where {I, ΩT <: Ω{I}, TAGS} =
    new{I, TAGS, ΩT}(ω, tags)
end

hastags(::Type{TaggedΩ{I, TAGS, ΩT}}, tags...) where {I, TAGS, ΩT} = hastags(TAGS, tags...)

tag(ω, tag_::Tags) = TaggedΩ(ω, tag_)
tag(tω::TaggedΩ, tag_::Tags) = TaggedΩ(tω.taggedω, merge(combinetag, tag_, tω.tags))

# Pass-throughs (tω::TaggedΩ should work like its tω.taggedω, but preserve tags)

Base.values(tω::TaggedΩ) = values(tω.taggedω)
Base.keys(tω::TaggedΩ) = keys(tω.taggedω)

@inline Base.getindex(tω::TaggedΩ, i) = TaggedΩ(getindex(tω.taggedω, i), tω.tags)
# Base.rand(tω::TaggedΩ, args...) = rand(tω.taggedω, args...; rng = rng(tω))

# Just pass through rand to tω.taggedω (have to repeat due to ambiguities with Random)
@inline Base.rand(tω::TaggedΩ) = rand(tω.taggedω; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, dims::Dims) = rand(tω.taggedω, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, dim::Integer, dims::Integer...) = rand(tω.taggedω, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, ::Type{T}) where T = rand(tω.taggedω, T; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ::Type{T}, dims::Dims) where T = rand(tω.taggedω, T, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ::Type{T}, dim::Integer, dims::Integer...) where T= rand(tω.taggedω, T, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, arr::Array) = rand(tω.taggedω, arr; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, arr::Array, dims::Dims) = rand(tω.taggedω, arr, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, arr::Array, dim::Integer, dims::Integer...) = rand(tω.taggedω, arr, dim, dims...; rng = rng(tω))

@inline Base.rand(tω::TaggedΩ, ur::UnitRange) = rand(tω.taggedω, ur; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ur::UnitRange, dims::Dims) = rand(tω.taggedω, ur, dims; rng = rng(tω))
@inline Base.rand(tω::TaggedΩ, ur::UnitRange, dim::Integer, dims::Integer...) = rand(tω.taggedω, ur, dim, dims...; rng = rng(tω))

@inline parentω(tω::TaggedΩ) = TaggedΩ(parentω(tω), tω.tags)


text:
α = 0.025

@testset "transform marginals" begin
  x1 = [0.2 0.4; 0.4 0.6; 0.6 0.8]
  convertmarg!(x1, TDist, [[10],[10]])
  @test x1 ≈ [-0.879058  -0.260185; -0.260185 0.260185; 0.260185 0.879058] atol=1.0e-5
  Random.seed!(43)
  x = rand(10000, 2)
  convertmarg!(x, Normal, [[0., 2.],[0., 3.]])
  @test pvalue(ExactOneSampleKSTest(x[:,1],Normal(0,2))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2],Normal(0,3))) > α
  Random.seed!(43)
  @test_throws AssertionError convertmarg!(randn(1000, 2), Normal)
end

@testset "generate corelation matrix" begin
  Random.seed!(43)
  @test cormatgen(2) ≈ [1.0 0.660768; 0.660768 1.0] atol=1.0e-5

  @test cormatgen_constant(3, 0.3) == [1 0.3 0.3; 0.3 1 0.3; 0.3 0.3 1]
  @test cormatgen_toeplitz(3, 0.3) == [1 0.3 0.09; 0.3 1 0.3; 0.09 0.3 1]
  @test cormatgen_two_constant(4, 0.5, 0.2) == [1. .5 .2 .2; .5 1. .2 .2; .2 .2 1. .2; .2 .2 .2 1.]

  @test cormatgen_constant_noised(10, 0.2, ϵ=0.) == cormatgen_constant(10, 0.2)
  @test all(eigvals(cormatgen_constant_noised(10, 0.2)) .> 0)
  @test diag(cormatgen_constant_noised(10, 0.2)) ≈ ones(10)

  @test cormatgen_toeplitz_noised(10, 0.2, ϵ=0.) == cormatgen_toeplitz(10, 0.2)
  @test all(eigvals(cormatgen_toeplitz_noised(10, 0.2)) .> 0)
  @test diag(cormatgen_toeplitz_noised(10, 0.2)) ≈ ones(10)

  @test cormatgen_two_constant_noised(10, 0.8, 0.2, ϵ=0.) == cormatgen_two_constant(10, 0.8, 0.2)
  @test all(eigvals(cormatgen_two_constant_noised(10, 0.8, 0.2)) .> 0)
  @test diag(cormatgen_two_constant_noised(10, 0.8, 0.2)) ≈ ones(10)

  Random.seed!(43)
  @test cormatgen_rand(2) ≈ [1.0 0.879086; 0.879086 1.0] atol=1.0e-5
  @test all(eigvals(cormatgen_rand(10)) .> 0)
  @test diag(cormatgen_rand(10)) ≈ ones(10)

end

@testset "gaussian copula" begin
  Random.seed!(43)
  x = gausscopulagen(500000, [1. 0.5; 0.5 1.])
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test tail(x[:,1], x[:,2], "l", 0.00001) ≈ 0
  @test tail(x[:,1], x[:,2], "r", 0.00001) ≈ 0
end
@testset "t-student copula" begin
  ν = 10
  rho = 0.5
  λ = 2*pdf(TDist(ν+1), -sqrt.((ν+1)*(1-rho)/(1+rho)))
  Random.seed!(43)
  xt = tstudentcopulagen(500000, [1. rho; rho 1.], ν);
  @test pvalue(ExactOneSampleKSTest(xt[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(xt[:,2], Uniform(0,1))) > α
  @test tail(xt[:,1], xt[:,2], "l") ≈ λ atol=1.0e-1
  @test tail(xt[:,1], xt[:,2], "r") ≈ λ atol=1.0e-1
  convertmarg!(xt, Normal)
  @test cor(xt) ≈ [1. rho; rho 1.] atol=1.0e-2
  Random.seed!(43)
  xtt = tstudentcopulagen(500000, [4. 2*rho; 2*rho 1.], ν);
  @test pvalue(ExactOneSampleKSTest(xtt[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(xtt[:,2], Uniform(0,1))) > α
  @test tail(xtt[:,1], xtt[:,2], "l") ≈ λ atol=1.0e-1
  @test tail(xtt[:,1], xtt[:,2], "r") ≈ λ atol=1.0e-1
  convertmarg!(xtt, Normal)
  @test cov(xtt) ≈ [1. rho; rho 1.] atol=1.0e-2
end
@testset "Frechet copula" begin
  Random.seed!(43)
  x = frechetcopulagen(500000, 3, 0.3);
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
  @test tail(x[:,1], x[:,2], "l") ≈ 0.3 atol=1.0e-1
  @test tail(x[:,1], x[:,2], "r") ≈ 0.3 atol=1.0e-1
  @test corspearman(x) ≈ [1. 0.3 0.3; 0.3 1. 0.3; 0.3 0.3 1.] atol=1.0e-2
  x = frechetcopulagen(500000, 3, 1.)
  ret = (x[:, 1] .<  0.2).*(x[:, 2] .<  0.3).*(x[:, 3] .<  0.5)
  @test length(findall(ret))/size(x,1) ≈ 0.2 atol=1.0e-3
  ret = (x[:, 1] .<  0.8).*(x[:, 2] .<  0.3).*(x[:, 3] .<  0.5)
  @test length(findall(ret))/size(x,1) ≈ 0.3 atol=1.0e-3
  ret = (x[:, 1] .<  0.8).*(x[:, 2] .<  0.9).*(x[:, 3] .<  0.5)
  @test length(findall(ret))/size(x,1) ≈ 0.5 atol=1.0e-3
  Random.seed!(43)
  x = frechetcopulagen(500000, 2, .0)
  @test corspearman(x) ≈ [1. 0.; 0. 1.] atol=1.0e-2
  x = frechetcopulagen(500000, 2, 0.8, 0.1);
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test corspearman(x) ≈ [1. 0.7; 0.7 1.] atol=1.0e-2
  @test tail(x[:,1], x[:,2], "l") ≈ 0.8 atol=1.0e-1
  @test tail(x[:,1], x[:,2], "r") ≈ 0.8 atol=1.0e-1
end

@testset "Frechet copula exceptions" begin
  @test_throws AssertionError frechetcopulagen(1000, 3, 0.8, 0.1)
  @test_throws DomainError frechetcopulagen(1000, 2, 0.8, 0.5)
  @test_throws DomainError frechetcopulagen(1000, 3, -0.8)
end

@testset "Marshall-Olkin helpers" begin
  @test τ2λ([0.4], [4.2, 1.2]) ≈ [4.2, 1.2, 3.6]
  @test τ2λ([0.5, 0.6, 0.5], [1., 2., 3., 4.]) ≈ [1, 2, 3, 0, 0.5, 0, 4]
  @test moρ2τ(0.6) ≈ 0.5 atol=1.0e-2
end
@testset "Marshall-Olkin copula" begin
  m = [0.252982 0.201189;  0.464758 0.409039; 0.585662 0.5357]
  @test mocopula([0.2 0.3 0.4; 0.3 0.4 0.6; 0.4 0.5 0.7], 2, [1., 1.5, 2.]) ≈ m atol=1.0e-4
  Random.seed!(43)
  x = marshallolkincopulagen(100000, [1.1, 0.2, 0.6])
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  a2 = 0.6/0.8
  a1 = 0.6/1.7
  @test corkendall(x)[1,2]≈ a1*a2/(a1+a2-a1*a2) atol=1.0e-3
  @test tail(x[:,1], x[:,2], "r") ≈ a1 atol=1.0e-1
  Random.seed!(42)
  x = marshallolkincopulagen(100000, [1.1, 0.2, 2.1, 0.6, 0.5, 3.2, 7.1, 2.1])
  @test pvalue(ExactOneSampleKSTest(x[:,1], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,2], Uniform(0,1))) > α
  @test pvalue(ExactOneSampleKSTest(x[:,3], Uniform(0,1))) > α
end