text:
group = "all"

if group == "fields" || group == "all"
    println(1)
    include("test_fields.jl")
end

if group == "grids" || group == "all"
    include("test_grids.jl")
end

if group == "melt" || group == "all"
    include("test_melt.jl")
end

if group == "models" || group == "all"
    include("test_models.jl")
end

if group == "outputs" || group == "all"
    include("test_outputting.jl")
end

if group == "simulations" || group == "all"
    include("test_simulations.jl")
end

if group == "timesteppingparams" || group == "all"
    include("test_timesteppingparams.jl")
end

if group == "utils" || group == "all"
    include("test_utils.jl")
end

if group == "verification" || group == "all"
    include("verification_tests.jl")
end

if group == "version_update" || group == "all"
    include(joinpath("version_update_test_verification","test_version_updates.jl"))
end


text:
"""
    y = bunquantize(q)

Calculate the value corresponding to a bidirectionally quantized quantity.
`q` is a 2n by m matrix containing the powers of 2 and their signs for each
quantized value.
"""
function bunquantize(q)
    n = size(q, 1) ÷ 2
    y = zeros(n)
    signs = 2:2:2*n
    powers = signs .- 1
    for i = 1:size(q, 2)
        y = y + 2.0 .^ q[powers,i] .* q[signs,i]
    end
    return y
end

ext: jl

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
#using HydrationOfCement
using Test

@testset "HydrationOfCement.jl" begin
    @test HydrationOfCement.coolFunc(2) == 4

end


text:
function tcm_predictHALOMNL(N, u, alpha_val, assortments)
	M = size(assortments,1)
	v = zeros(size(assortments))
	for m = 1 : M
		v[m,N] = 1.0
		for p in 1:N-1
			if (assortments[m,p] > 0)
				v[m,p] = exp(u[p] + dot( (1 - assortments[m,1:N-1]), alpha_val[:, p] )   )
			end
		end
		v[m,:] = v[m,:] / sum(v[m,:])
	end

	return v
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
# Spatial calculations

export
    vector1D,
    vector,
    wrap_coords,
    wrap_coords_vec,
    maxwell_boltzmann,
    random_velocities,
    random_velocities!,
    bond_angle,
    torsion_angle

"""
    vector1D(c1, c2, side_length)

Displacement between two 1D coordinate values from c1 to c2, accounting for
the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinate accounting for the periodic boundaries.
"""
function vector1D(c1, c2, side_length)
    if c1 < c2
        return (c2 - c1) < (c1 - c2 + side_length) ? (c2 - c1) : (c2 - c1 - side_length)
    else
        return (c1 - c2) < (c2 - c1 + side_length) ? (c2 - c1) : (c2 - c1 + side_length)
    end
end

"""
    vector(c1, c2, box_size)

Displacement between two coordinate values from c1 to c2, accounting for
the bounding box.
The minimum image convention is used, so the displacement is to the closest
version of the coordinates accounting for the periodic boundaries.
"""
vector(c1, c2, box_size) = vector1D.(c1, c2, box_size)

@generated function vector(c1::SVector{N}, c2::SVector{N}, box_size) where N
    quote
        Base.Cartesian.@ncall $N SVector{$N} i->vector1D(c1[i], c2[i], box_size[i])
    end
end

square_distance(i, j, coords, box_size) = sum(abs2, vector(coords[i], coords[j], box_size))

# Pad a vector to 3D to allow operations such as the cross product
function vector_pad3D(c1::SVector{2, T}, c2::SVector{2, T}, box_size::SVector{2, T}) where T
    SVector{3, T}(
        vector1D(c1[1], c2[1], box_size[1]),
        vector1D(c1[2], c2[2], box_size[2]),
        zero(T),
    )
end

vector_pad3D(c1::SVector{3}, c2::SVector{3}, box_size::SVector{3}) = vector(c1, c2, box_size)

# Trim a vector back to 2D if required
trim3D(v::SVector{3, T}, box_size::SVector{2}) where T = SVector{2, T}(v[1], v[2])
trim3D(v::SVector{3}, box_size::SVector{3}) = v

"""
    wrap_coords(c, side_length)

Ensure a 1D coordinate is within the simulation box and return the coordinate.
"""
wrap_coords(c, side_length) = c - floor(c / side_length) * side_length

"""
    wrap_coords_vec(c, box_size)

Ensure a coordinate is within the simulation box and return the coordinate.
"""
wrap_coords_vec(v, box_size) = wrap_coords.(v, box_size)

"""
    velocity(mass, temperature; dims=3)

Generate a random velocity from the Maxwell-Boltzmann distribution.
"""
function AtomsBase.velocity(mass, temp; dims::Integer=3)
    return SVector([maxwell_boltzmann(mass, temp) for i in 1:dims]...)
end

const mb_conversion_factor = uconvert(u"u * nm^2 * ps^-2 * K^-1", Unitful.k)

"""
    maxwell_boltzmann(mass, temperature)

Generate a random speed along one dimension from the Maxwell-Boltzmann distribution.
"""
function maxwell_boltzmann(mass, temp)
    T = typeof(convert(AbstractFloat, ustrip(temp)))
    k = unit(temp) == NoUnits ? T(ustrip(mb_conversion_factor)) : T(mb_conversion_factor)
    σ = sqrt(k * temp / mass)
    return rand(Normal(zero(T), T(ustrip(σ)))) * unit(σ)
end

"""
    random_velocities(sys, temp)

Generate random velocities from the Maxwell-Boltzmann distribution
for a `System`.
"""
function random_velocities(sys::System{D}, temp) where D
    [velocity(a.mass, temp; dims=D) for a in Array(sys.atoms)]
end

"""
    random_velocities!(sys, temp)

Set the velocities of a `System` to random velocities generated from the
Maxwell-Boltzmann distribution.
"""
function random_velocities!(sys::System, temp)
    sys.velocities = random_velocities(sys, temp)
    return sys
end

# Sometimes domain error occurs for acos if the value is > 1.0 or < -1.0
acosbound(x::Real) = acos(clamp(x, -1, 1))

"""
    bond_angle(coord_i, coord_j, coord_k, box_size)
    bond_angle(vec_ba, vec_bc)

Calculate the bond or pseudo-bond angle in radians between three
coordinates or two vectors.
The angle between B→A and B→C is returned in the range 0 to π.
"""
function bond_angle(coords_i, coords_j, coords_k, box_size)
    vec_ba = vector(coords_j, coords_i, box_size)
    vec_bc = vector(coords_j, coords_k, box_size)
    return bond_angle(vec_ba, vec_bc)
end

function bond_angle(vec_ba, vec_bc)
    acosbound(dot(vec_ba, vec_bc) / (norm(vec_ba) * norm(vec_bc)))
end

"""
    torsion_angle(coord_i, coord_j, coord_k, coord_l, box_size)
    torsion_angle(vec_ab, vec_bc, vec_cd)

Calculate the torsion angle in radians defined by four coordinates or
three vectors.
The angle between the planes defined by atoms (i, j, k) and (j, k, l) is
returned in the range -π to π.
"""
function torsion_angle(coords_i, coords_j, coords_k, coords_l, box_size)
    vec_ab = vector(coords_i, coords_j, box_size)
    vec_bc = vector(coords_j, coords_k, box_size)
    vec_cd = vector(coords_k, coords_l, box_size)
    return torsion_angle(vec_ab, vec_bc, vec_cd)
end

function torsion_angle(vec_ab, vec_bc, vec_cd)
    cross_ab_bc = vec_ab × vec_bc
    cross_bc_cd = vec_bc × vec_cd
    θ = atan(
        ustrip(dot(cross_ab_bc × cross_bc_cd, normalize(vec_bc))),
        ustrip(dot(cross_ab_bc, cross_bc_cd)),
    )
    return θ
end


text:
# NOTE: We want to load only (part of) one momentum sector of the diagonalized
# Hamiltonian into memory at a time.  We therefore make two passes over all the
# momenta: once when constructing the original state in the energy basis, and
# again when moving the time evolved states back to the position basis.
#
# (FIXME: Actually, it is almost certainly possible to do time evolution using
# just a single pass over the momentum eigenstates, though this may involve
# giving up some flexibility.  Also, for systems that have time-reversal
# invariant and invariant under k <-> -k, we should be able to load eigenstates
# corresponding to only half the Brillouin zone.)

# Implementation of matrix multiplication for JldDataset.  We load one column
# at a time since the matrix is stored in column-major form.  See discussion at
# https://groups.google.com/d/msg/julia-users/cEHmqsMLX5c/NXREyAsVBwAJ

function my_Ac_mul_B(a::JLD.JldDataset, b::Matrix)
    K = size(b, 2)
    M = size(a, 2)
    N = size(a, 1)
    N == size(b, 1) || throw(DimensionMismatch())
    rv = zeros(eltype(b), M, K) # XXX: would be better to use promote_type. see https://github.com/JuliaLang/JLD.jl/issues/28
    for j in 1:M
        col = a[:,j]
        for k in 1:K
            @simd for i in 1:N
                # NOTE: col[i,1] is necessary due to
                # https://github.com/JuliaLang/HDF5.jl/issues/267
                @inbounds rv[j,k] += conj(col[i,1]) * b[i,k]
            end
        end
    end
    return rv
end

function my_A_mul_B!(c::Matrix, a::JLD.JldDataset, b::Matrix)
    @assert c !== b
    K = size(b, 2)
    M = size(a, 1)
    N = size(a, 2)
    size(b, 1) == N || throw(DimensionMismatch())
    size(c) == (M, K) || throw(DimensionMismatch())
    fill!(c, zero(eltype(c)))
    for j in 1:N
        col = a[:,j]
        for k in 1:K
            @simd for i in 1:M
                # NOTE: col[i,1] is necessary due to
                # https://github.com/JuliaLang/HDF5.jl/issues/267
                @inbounds c[i,k] += b[j,k] * col[i,1]
            end
        end
    end
    return c
end

# These functions should behave as usual for all other data types (e.g. normal
# vectors/matrices)
my_Ac_mul_B(a, b) = a' * b
my_A_mul_B!(c, a, b) = mul!(c, a, b)

function calculate_energy_basis_size(state_table::RepresentativeStateTable, k_indices)
    if k_indices == eachmomentumindex(state_table.hs.lattice)
        return length(state_table.hs.indexer)
    end

    energy_basis_size = 0
    for sector_index in 1:state_table.sector_count
        for momentum_index in k_indices
            energy_basis_size += length(DiagonalizationSector(state_table, sector_index, momentum_index))
        end
    end
    return energy_basis_size
end

function to_energy_basis(load_momentum_sector::Function, state_table::RepresentativeStateTable, initial_states::VecOrMat; k_indices=eachmomentumindex(state_table.hs.lattice))
    if size(initial_states, 1) != length(state_table.hs.indexer)
        throw(ArgumentError("Initial state must match indexer size"))
    end
    nstates = size(initial_states, 2)

    energy_basis_size = calculate_energy_basis_size(state_table, k_indices)

    ###
    # Transform initial state to energy basis
    ###
    initial_energy_states = Array{ComplexF64}(undef, energy_basis_size, size(initial_states)[2:end]...)
    all_energies = sizehint!(Float64[], energy_basis_size)
    offset = 0
    initial_momentum_state = ComplexF64[]
    for sector_index in 1:state_table.sector_count
        for momentum_index in k_indices
            load_momentum_sector(sector_index, momentum_index) do reduced_indexer, reduced_energies, reduced_eigenstates
                @assert length(reduced_indexer) == length(reduced_energies) == size(reduced_eigenstates, 1) == size(reduced_eigenstates, 2)
                myrange = offset+1 : offset+length(reduced_indexer)
                diagsect = DiagonalizationSector(state_table, sector_index, momentum_index, reduced_indexer)

                initial_momentum_states = zeros(ComplexF64, length(diagsect), nstates)

                # Project each state onto current momentum basis
                #
                # FIXME: actually, i should make a function in abelian.jl that does
                # this (currently apply_reduced_hamiltonian is the only similar
                # thing).  if so, make sure it fills with zeros first.
                for x in 1:nstates
                    for (i, (reduced_i, alpha)) in enumerate(diagsect.representative_v)
                        if reduced_i != 0
                            initial_momentum_states[reduced_i, x] += initial_states[i, x] * conj(alpha)
                        end
                    end
                end

                # Transform to energy basis
                initial_energy_states[myrange,:] = my_Ac_mul_B(reduced_eigenstates, initial_momentum_states)

                append!(all_energies, reduced_energies)

                offset += length(reduced_indexer)
                @assert length(all_energies) == offset
                nothing
            end
        end
    end
    @assert offset == energy_basis_size

    return initial_energy_states, all_energies
end

function time_evolve_to_position_basis(load_momentum_sector::Function, state_table::RepresentativeStateTable, initial_energy_states::VecOrMat, time_steps::AbstractVector{<:Real}; k_indices=eachmomentumindex(state_table.hs.lattice))
    energy_basis_size = calculate_energy_basis_size(state_table, k_indices)
    if size(initial_energy_states, 1) != energy_basis_size
        throw(ArgumentError("Initial energy state must match indexer size"))
    end

    ###
    # Time evolve and move back to position basis
    ###
    basis_size = length(state_table.hs.indexer)
    output_states = zeros(ComplexF64, basis_size, length(time_steps), size(initial_energy_states)[2:end]...)
    offset = 0
    for sector_index in 1:state_table.sector_count
        for momentum_index in k_indices
            load_momentum_sector(sector_index, momentum_index) do reduced_indexer, reduced_energies, reduced_eigenstates
                @assert length(reduced_indexer) == length(reduced_energies) == size(reduced_eigenstates, 1) == size(reduced_eigenstates, 2)
                diagsect = DiagonalizationSector(state_table, sector_index, momentum_index, reduced_indexer)

                momentum_states = Array{ComplexF64}(undef, length(diagsect), length(time_steps))
                time_evolved_sector = Array{ComplexF64}(undef, length(diagsect), length(time_steps))

                # Loop through each initial state
                for z in 1:size(initial_energy_states, 2)
                    for (t_i, t) in enumerate(time_steps)
                        # Time evolve
                        for i in 1:length(reduced_indexer)
                            time_evolved_sector[i, t_i] = initial_energy_states[offset + i, z] * exp(-im * reduced_energies[i] * t)
                        end
                    end

                    # Move back to momentum basis
                    #
                    # FIXME: currently this means loading the reduced_eigenstates
                    # once for each initial state.  It would be great to be able to
                    # do it once entirely.
                    my_A_mul_B!(momentum_states, reduced_eigenstates, time_evolved_sector)

                    for (t_i, t) in enumerate(time_steps)
                        # Move back to position basis
                        for (reduced_i, i, alpha) in diagsect.coefficient_v
                            output_states[i, t_i, z] += momentum_states[reduced_i, t_i] * alpha
                        end
                    end
                end

                offset += length(reduced_indexer)
                nothing
            end
        end
    end
    @assert offset == energy_basis_size

    return output_states
end

function time_evolve(load_momentum_sector::Function, state_table::RepresentativeStateTable, initial_states::VecOrMat, time_steps::AbstractVector{Float64}; kwargs...)
    ψ_e, = to_energy_basis(load_momentum_sector, state_table, initial_states; kwargs...)
    return time_evolve_to_position_basis(load_momentum_sector, state_table, ψ_e, time_steps; kwargs...)
end


text:
# # [Quantum Fourier Transformation and Phase Estimation](@id example-qft)

# Let's use Yao first

using Yao

# ## Quantum Fourier Transformation

# The Quantum Fourier Transformation (QFT) circuit is to repeat
# two kinds of blocks repeatly:

# ![qft-circuit](assets/qft.png)

# The basic building block control phase shift gate is defined
# as

# ```math
# R(k)=\begin{bmatrix}
# 1 & 0\\
# 0 & \exp\left(\frac{2\pi i}{2^k}\right)
# \end{bmatrix}
# ```

# Let's define block `A` and block `B`, block `A` is actually
# a control block.

A(i, j) = control(i, j=>shift(2π/(1<<(i-j+1))))

# Once you construct the blockl you can inspect its matrix using `mat`
# function. Let's construct the circuit in dash box A, and see the matrix of
# ``R_4`` gate.

R4 = A(4, 1)

# If you have read about [preparing GHZ state](@ref example-ghz),
# you probably know that in Yao, we could just leave the number of qubits, and it
# will be evaluated when possible.

R4(5)

# its matrix will be

mat(R4(5))

# Then we repeat this control block over
# and over on different qubits, and put a Hadamard gate
# to `i`th qubit to construct `i`-th `B` block.

B(n, k) = chain(n, j==k ? put(k=>H) : A(j, k) for j in k:n)

# We need to input the total number of qubits `n` here because we have to iterate
# through from `k`-th location to the last.

# Now, let's construct the circuit by chaining all the `B` blocks together

qft(n) = chain(B(n, k) for k in 1:n)
qft(4)

# ## Wrap QFT to an external block

# In most cases, `function`s are enough to wrap quantum circuits, like `A`
# and `B` we defined above, but sometimes, we need to dispatch specialized
# methods on certain kinds of quantum circuit, or we want to define an external
# block to export, thus, it's useful to be able to wrap circuit to custom blocks.

# First, we define a new type as subtype of `PrimitiveBlock` since we are not
# going to use the subblocks of `QFT`, if you need to use its subblocks, it'd
# be better to define it under `CompositeBlock`.

struct QFT{N} <: PrimitiveBlock{N} end
QFT(n::Int) = QFT{n}()

# Now, let's define its circuit

circuit(::QFT{N}) where N = qft(N)

# And forward `mat` to its circuit's matrix

YaoBlocks.mat(::Type{T}, x::QFT) where T = mat(T, circuit(x))

# You may notice, it is a little ugly to print `QFT` at the moment,
# this is because we print the type summary by default, you can define
# your own printing by overloading `print_block`

YaoBlocks.print_block(io::IO, x::QFT{N}) where N = print(io, "QFT($N)")

# Since it is possible to use FFT to simulate the results of QFT (like cheating),
# we could define our custom `apply!` method:

using FFTW, LinearAlgebra

function YaoBlocks.apply!(r::ArrayReg, x::QFT)
    α = sqrt(length(statevec(r)))
    invorder!(r)
    lmul!(α, ifft!(statevec(r)))
    return r
end

# Now let's check if our `apply!` method is correct:

r = rand_state(5)
r1 = r |> copy |> QFT(5)
r2 = r |> copy |> circuit(QFT(5))
r1 ≈ r2

# We can get iQFT (inverse QFT) directly by calling `adjoint`

QFT(5)'

# QFT and iQFT are different from FFT and IFFT in three ways,

# 1. they are different by a factor of ``\sqrt{2^n}`` with ``n`` the number of qubits.
# 2. the [bit numbering](https://quantumbfs.github.io/BitBasis.jl/stable/tutorial/#Conventions-1) will exchange after applying QFT or iQFT.
# 3. due to the convention, QFT is more related to IFFT rather than FFT.

# ## Phase Estimation

# Since we have QFT and iQFT blocks we can then use them to
# realize phase estimation circuit, what we want to realize
# is the following circuit:

# ![phase estimation](assets/phaseest.png)

# First we call Hadamard gates repeatly on first `n` qubits.

Hadamards(n) = repeat(H, 1:n)

# Then in dashed box `B`, we have controlled unitaries:

ControlU(n, m, U) = chain(n+m, control(k, n+1:n+m=>matblock(U^(2^(k-1)))) for k in 1:n)

# each of them is a `U` of power ``2^{(k-1)}``.

# Since we will only apply the qft and Hadamard on first `n` qubits,
# we could use `Subroutine`, which creates a context of
# a sub-scope of the qubits.

PE(n, m, U) =
    chain(n+m, # total number of the qubits
        subroutine(Hadamards(n), 1:n), # apply H in local scope
        ControlU(n, m, U),
        subroutine(QFT(n)', 1:n))

# we use the first `n` qubits as the output space to store phase ``ϕ``, and the
# other `m` qubits as the input state which corresponds to an eigenvector of
# oracle matrix `U`.

# The subroutine here uses `focus!` and `relax!` to manage
# a local scope of quantum circuit, and only activate the first `n` qubits while applying
# the block inside the subroutine context, and the scope will be `relax!`ed
# back, after the context. This is equivalent to manually `focus!`
# then `relax!`

# fullly activated

r = rand_state(5)

# first 3 qubits activated

focus!(r, 1:3)

# relax back to the original

relax!(r, 1:3)

# In this way, we will be able to apply small operator directly
# on the subset of the qubits.

# Details about the algorithm can be found here:
# [Quantum Phase Estimation Algorithm](https://en.wikipedia.org/wiki/Quantum_phase_estimation_algorithm)

# Now let's check the results of our phase estimation.

# First we need to set up a unitary with known phase, we set the phase to be
# 0.75, which is `0.75 * 2^3 == 6 == 0b110` .

N, M = 3, 5
P = eigen(rand_unitary(1<<M)).vectors
θ = Int(0b110) / 1<<N
phases = rand(1<<M)
phases[0b010+1] = θ
U = P * Diagonal(exp.(2π * im * phases)) * P'

# and then generate the state ``ψ``

psi = P[:, 3]

# In the phase estimation process, we will feed the state to circuit and measure
# the first `n` qubits processed by iQFT.

r = join(ArrayReg(psi), zero_state(N))
r |> PE(N, M, U)

# Since our phase can be represented by 3 qubits precisely, we only need to measure once

results = measure(r, 1:N; nshots=1)

# Recall that our QFT's bit numbering is reversed, let's reverse it back

using BitBasis
estimated_phase = bfloat(results[]; nbits=N)

# the phase is exactly `0.75`!


text:
struct DirectSummand{T, M<:AbstractMatrix{T}} <: AbstractMatrix{T}
    basis::M
    multiplicity::Int
    degree::Int
    simple::Bool
end

image_basis(ds::DirectSummand) = ds.basis
issimple(ds::DirectSummand) = ds.simple
PermutationGroups.degree(ds::DirectSummand) = ds.degree
multiplicity(ds::DirectSummand) = ds.multiplicity

Base.size(ds::DirectSummand) = size(ds.basis)
Base.@propagate_inbounds Base.getindex(ds::DirectSummand, i...) =
    image_basis(ds)[i...]

function SparseArrays.sparse(ds::DirectSummand)
    sp = sparse(image_basis(ds))
    return DirectSummand(sp, multiplicity(ds), degree(ds), issimple(ds))
end

function SparseArrays.droptol!(
    ds::DirectSummand{T, <:AbstractSparseArray},
    tol
) where T
    droptol!(image_basis(ds), tol)
    return ds
end


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

struct CustomMetricsV1beta1Api <: SwaggerApi
    client::Swagger.Client
end

"""

Retrieve the given metric for the given non-namespaced object (e.g. Node, PersistentVolume). Composite metric name is of the form \"<objecttype>/<objectname>/<metricname>\" Passing '*' for objectname would retrieve the given metric for all non-namespaced objects of the given type.
Param: compositemetricname::String (required)
Param: pretty::String
Param: fieldSelector::String
Param: labelSelector::String
Param: limit::Int32
Return: IoK8sApiCustomMetricsV1beta1MetricValueList
"""
function listCustomMetricsV1beta1MetricValue(_api::CustomMetricsV1beta1Api, compositemetricname::String; pretty=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApiCustomMetricsV1beta1MetricValueList, "/apis/custom.metrics.k8s.io/v1beta1/{compositemetricname}", ["BearerToken"])
    Swagger.set_param(_ctx.path, "compositemetricname", compositemetricname)  # type String
    Swagger.set_param(_ctx.query, "pretty", pretty)  # type String
    Swagger.set_param(_ctx.query, "fieldSelector", fieldSelector)  # type String
    Swagger.set_param(_ctx.query, "labelSelector", labelSelector)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf", "application/json;stream=watch", "application/vnd.kubernetes.protobuf;stream=watch"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

retrieve the given metric (in composite form) which describes the given namespace. Composite form of metrics can be either \"metric/<metricname>\" to fetch metrics of all objects in the namespace, or \"<objecttype>/<objectname>/<metricname>\" to fetch metrics of the specified object type and name. Passing \"*\" for objectname would fetch metrics of all objects of the specified type in the namespace.
Param: compositemetricname::String (required)
Param: namespace::String (required)
Param: pretty::String
Param: fieldSelector::String
Param: labelSelector::String
Param: limit::Int32
Return: IoK8sApiCustomMetricsV1beta1MetricValueList
"""
function listCustomMetricsV1beta1NamespacedMetricValue(_api::CustomMetricsV1beta1Api, compositemetricname::String, namespace::String; pretty=nothing, fieldSelector=nothing, labelSelector=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IoK8sApiCustomMetricsV1beta1MetricValueList, "/apis/custom.metrics.k8s.io/v1beta1/namespaces/{namespace}/{compositemetricname}", ["BearerToken"])
    Swagger.set_param(_ctx.path, "compositemetricname", compositemetricname)  # type String
    Swagger.set_param(_ctx.path, "namespace", namespace)  # type String
    Swagger.set_param(_ctx.query, "pretty", pretty)  # type String
    Swagger.set_param(_ctx.query, "fieldSelector", fieldSelector)  # type String
    Swagger.set_param(_ctx.query, "labelSelector", labelSelector)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json", "application/yaml", "application/vnd.kubernetes.protobuf", "application/json;stream=watch", "application/vnd.kubernetes.protobuf;stream=watch"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["*/*"] : [_mediaType])
    Swagger.exec(_ctx)
end

export listCustomMetricsV1beta1MetricValue, listCustomMetricsV1beta1NamespacedMetricValue


text:
using StatsBase
using Base.Test

# test hist
@test sum(fit(Histogram,[1,2,3]).weights) == 3
@test fit(Histogram,Int[]).weights == Int[]
@test fit(Histogram,[1]).weights == [1]
@test fit(Histogram,[1,2,3],[0,2,4]) == Histogram([0,2,4],[2,1])
@test fit(Histogram,[1,2,3],[0,2,4]) != Histogram([0,2,4],[1,1])
@test fit(Histogram,[1,2,3],0:2:4) == Histogram(0:2:4,[2,1])
@test all(fit(Histogram,[1:100;]/100,0.0:0.01:1.0).weights .==1)
@test fit(Histogram,[1,1,1,1,1]).weights[1] == 5
@test sum(fit(Histogram,(rand(100),rand(100))).weights) == 100
@test fit(Histogram,1:100,nbins=5,closed=:right).weights == [20,20,20,20,20]
@test fit(Histogram,1:100,nbins=5,closed=:left).weights == [19,20,20,20,20,1]
@test fit(Histogram,0:99,nbins=5,closed=:right).weights == [1,20,20,20,20,19]
@test fit(Histogram,0:99,nbins=5,closed=:left).weights == [20,20,20,20,20]

@test fit(Histogram,(1:100,1:100),nbins=5).weights == diagm([20,20,20,20,20])
@test fit(Histogram,(1:100,1:100),nbins=(5,5)).weights == diagm([20,20,20,20,20])

@test fit(Histogram,1:100,weights(ones(100)),nbins=5).weights == [20,20,20,20,20]
@test fit(Histogram,1:100,weights(2*ones(100)),nbins=5).weights == [40,40,40,40,40]

@test eltype(fit(Histogram,1:100,weights(ones(Int,100)),nbins=5).weights) == Int
@test eltype(fit(Histogram,1:100,weights(ones(Float64,100)),nbins=5).weights) == Float64

# histrange
# Note: atm histrange must be qualified
@test StatsBase.histrange(Float64[], 0, :left) == 0.0:1.0:0.0
@test StatsBase.histrange(Float64[1:5;], 1, :left) == 0.0:5.0:10.0
@test StatsBase.histrange(Float64[1:10;], 1, :left) == 0.0:10.0:20.0
@test StatsBase.histrange(1.0, 10.0, 1, :left) == 0.0:10.0:20.0

@test StatsBase.histrange([0.201,0.299], 10, :left) == 0.2:0.01:0.3
@test StatsBase.histrange([0.2,0.299], 10, :left) == 0.2:0.01:0.3
@test StatsBase.histrange([0.2,0.3], 10, :left)  == 0.2:0.01:0.31
@test StatsBase.histrange(0.2, 0.3,  10, :left)  == 0.2:0.01:0.31
@test StatsBase.histrange([0.2,0.3], 10, :right) == 0.19:0.01:0.3
@test StatsBase.histrange(0.2, 0.3,  10, :right) == 0.19:0.01:0.3

@test StatsBase.histrange([200.1,299.9], 10, :left) == 200.0:10.0:300.0
@test StatsBase.histrange([200.0,299.9], 10, :left) == 200.0:10.0:300.0
@test StatsBase.histrange([200.0,300.0], 10, :left) == 200.0:10.0:310.0
@test StatsBase.histrange([200.0,300.0], 10, :right) == 190.0:10.0:300.0

@test StatsBase.histrange(Int64[1:5;], 1, :left) == 0:5:10
@test StatsBase.histrange(Int64[1:10;], 1, :left) == 0:10:20

@test StatsBase.histrange([1, 2, 3, 4], 4) == 0.0:1.0:4.0
@test StatsBase.histrange([1, 2, 2, 4], 4) == 0.0:1.0:4.0
@test StatsBase.histrange([1, 10], 4) == 0.0:5.0:10.0
@test StatsBase.histrange([1, 20], 4) == 0.0:5.0:20.0
@test StatsBase.histrange([1, 600], 4) == 0.0:200.0:600.0
@test StatsBase.histrange([1, -1000], 4) == -1500.0:500.0:500.0

# Base issue #13326
l,h = extrema(StatsBase.histrange([typemin(Int),typemax(Int)], 10))
@test l <= typemin(Int)
@test h >= typemax(Int)

@test_throws ArgumentError StatsBase.histrange([1, 10], 0)
@test_throws ArgumentError StatsBase.histrange([1, 10], -1)
@test_throws ArgumentError StatsBase.histrange([1.0, 10.0], 0)
@test_throws ArgumentError StatsBase.histrange([1.0, 10.0], -1)
@test_throws ArgumentError StatsBase.histrange(Float64[],-1)
@test_throws ArgumentError StatsBase.histrange([0.], 0)

# hist show
show_h = sprint(show, fit(Histogram,[1,2,3]))
@test contains(show_h, "edges:\n  0.0:1.0:3.0")
@test contains(show_h, "weights: $([1,1,1])")
@test contains(show_h, "closed: right")


text:
using Catalyst
using Distributions, Random
using DelaySSAToolkit

rn = @reaction_network begin
    A, 0 --> X
    B, Y --> 0
end A B

u0 = [0,0]
de_chan0 = [[]]
tf = 100.
tspan = (0,tf)
ps = [10., 1.] # dummy variables, later on will be drawn from a Gamma distribution
τ = 1 # dummy variables, later on will be drawn from a Gamma distribution
delay_trigger = Dict(1=>[1=>τ])
delay_complete = Dict(1=>[1=>-1, 2=>1])
delay_interrupt = Dict()

jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)
dprob = DiscreteProblem(jumpsys,u0,tspan,ps)
# use ensemble problem
algo = DelayRejection()
delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
djprob = DelayJumpProblem(jumpsys, dprob, algo, delayjumpset, de_chan0, save_positions=(false,false))

function prob_func1(prob, i, repeat)
    rng = Random.seed!(i)
    A = rand(rng, Gamma(8,1/0.23))
    B = rand(rng, Gamma(9,1/625))
    τ  = rand(rng, Gamma(7,1))
    remake(prob, p = [A,B], delay_trigger=Dict(1=>[1=>τ]))
end
function prob_func2(prob, i, repeat)
    rng = Random.seed!(i)
    α = rand(rng, Gamma(63,1/9)) # the Gamma distribution uses shape α, scale θ, Gamma(α,θ). In the paper, Gamma distribution uses shape α and rate  β. One needs to set the inverse.
    β = rand(rng, Gamma(10.,1/10.))
    A = rand(rng, Gamma(8,1/0.23))
    B = rand(rng, Gamma(9,1/625))
    τ  = rand(rng, Gamma(α,1/β))
    remake(prob, p = [A,B], delay_trigger=Dict(1=>[1=>τ]))
end

prob1 = EnsembleProblem(djprob, prob_func = prob_func1)
@time ens1 = solve(prob1, SSAStepper(), EnsembleThreads(), trajectories = 40, saveat= 1.)
Y_evolution  = [ens1[i][2,:] for i in 1:40]

using DifferentialEquations.EnsembleAnalysis
using Plots;theme(:vibrant)
plot(Y_evolution, linealpha=0.4, legend=false , linewidth = 2, fmt=:svg, xlabel = "Time", ylabel="# of Y individuals")
savefig("docs/src/assets/heterogeneous_delay1.svg")

prob2 = EnsembleProblem(djprob, prob_func = prob_func2)
@time ens2 = solve(prob2, SSAStepper(), EnsembleThreads(), trajectories = 10^4)

last_slice = componentwise_vectors_timepoint(ens2, tf)
histogram(last_slice[1], bins=0:20:1000, normalize=:pdf, xlabel = "# of X individuals", ylabel = "Probability",fillalpha = 0.6, linecolor = :orange, legend = false)
# savefig("docs/src/assets/heterogeneous_delay2.svg")


text:
"""
    geterror()

Returns the last slurm error number and the error string.
"""
function geterror()
  err = api.slurm_get_errno()
  return err, unsafe_string(api.slurm_strerror(err))
end


text:
export velocities, positions, pitches, durations
export replace_notes, getnotnotes, modpositions
export MoreVelNote, randomnotes
using MIDI: AbstractNote

"""
    MoreVelNote

Enables velocities higher than 127 but other than that identical with [`Note`](@ref).

When [`Notes`](@ref) with `MoreVelNote` are written into a [`MIDITrack`](@ref),
any velocity higher than 127 is equaled to 127.
"""
mutable struct MoreVelNote <: AbstractNote
    pitch::UInt8
    velocity::UInt
    position::UInt
    duration::UInt
    channel::UInt8

    MoreVelNote(pitch, velocity, position, duration, channel = 0) =
        if channel > 0x7F
            error( "Channel must be less than 128" )
        else
            new(pitch, velocity, position, duration, channel)
        end
end
Notes(::Type{MoreVelNote}) = Notes{MoreVelNote}(Vector{MoreVelNote}[], 960)

"""
    replace_notes(oldtrack::MIDITrack, notes::Notes) -> newtrack
Create a *new* midi track that copies all `MIDIEvent`s that
are not NOTEON or NOTEOFF from the
`oldtrack` while it also adds all notes from `notes`.
"""
function replace_notes(oldtrack::MIDI.MIDITrack, notes::Notes)

    newtrack = MIDI.MIDITrack()

    # Get pedal events with relative time and absolute time
    other_events_abspos, other_events = getnotnotes(oldtrack)

    # Second track will contain the new `notes`
    MIDI.addnotes!(newtrack, notes)

    addevents!(newtrack, other_events_abspos, other_events)

    return newtrack
end

"""
    getnotnotes(track::MIDI.MIDITrack) -> (abs_pos, events)
Find all events in `track` that are not NOTEON or NOTEOFF.
Return the found events and their positions in absolute time (in ticks).

Each event can be added to another `MIDITrack` using
```julia
for ev in zip(abs_pos, events)
    MIDI.addevent!(newtrack, ev...)
end
```

See also [`replace_notes`](@ref).
"""
function getnotnotes(oldtrack::MIDI.MIDITrack)
    old_events = oldtrack.events
    # Get pedal events with relative time and absolute time
    other_events = MIDI.MIDIEvent[]
    other_events_abspos = Int[]
    abstime = 0
    for i in 1:length(old_events)
        abstime += old_events[i].dT
        if typeof(old_events[i]) == MIDI.MIDIEvent
            # skip NOTEON and NOTEOFF events of first channel
            if old_events[i].status != 0x80 && old_events[i].status != 0x90
                push!(other_events, old_events[i])
                push!(other_events_abspos, abstime)
            end
        end
    end

    return other_events_abspos, other_events
end


text:
using FastIntegers
using Base.Test

for T in (Int, UInt)
    f1,f2 = map(T, T <: Unsigned ? (1,3) : (-1,+1))
    for x in (f1*42, f2*42)

        @test fast_pos(x) === +x
        @test fast_neg(x) === -x
        @test fast_abs(x) === abs(x)

        for y in (f1*13, f2*13)

            @test fast_add(x, y) === x + y
            @test fast_sub(x, y) === x - y
            @test fast_mul(x, y) === x * y
            @test fast_div(x, y) === x ÷ y
            @test fast_rem(x, y) === x % y
            @test fast_fld(x, y) === fld(x, y)
            @test fast_mod(x, y) === mod(x, y)
            @test fast_cld(x, y) === cld(x, y)

        end
    end
end

@test fast_abs(typemin(Int)) === typemin(Int)


text:
module ReactiveMPMessageTest

using Test
using ReactiveMP

import Base.Iterators: repeated, product

@testset "Message" begin

    @testset "Constructor" begin
        data    = PointMass(1)

        @test getdata(Message(data, true, true))    === data
        @test is_clamped(Message(data, true, true)) === true
        @test is_initial(Message(data, true, true)) === true

        @test getdata(Message(data, true, false))    === data
        @test is_clamped(Message(data, true, false)) === true
        @test is_initial(Message(data, true, false)) === false

        @test getdata(Message(data, false, true))    === data
        @test is_clamped(Message(data, false, true)) === false
        @test is_initial(Message(data, false, true)) === true

        @test getdata(Message(data, false, false))    === data
        @test is_clamped(Message(data, false, false)) === false
        @test is_initial(Message(data, false, false)) === false
    end

    @testset "multiply_messages" begin
        dist1 = NormalMeanVariance(randn(), rand())
        dist2 = NormalMeanVariance(randn(), rand())

        @test getdata(Message(dist1, false, false) * Message(dist2, false, false))    == prod(ProdAnalytical(), dist1, dist2)
        @test getdata(Message(dist2, false, false) * Message(dist1, false, false))    == prod(ProdAnalytical(), dist2, dist1)

        for (left_is_initial, right_is_initial) in product(repeated([ true, false ], 2)...)
            @test is_clamped(Message(dist1, true, left_is_initial) * Message(dist2, false, right_is_initial)) == false
            @test is_clamped(Message(dist1, false, left_is_initial) * Message(dist2, true, right_is_initial)) == false
            @test is_clamped(Message(dist1, true, left_is_initial) * Message(dist2, true, right_is_initial))  == true
            @test is_clamped(Message(dist2, true, left_is_initial) * Message(dist1, false, right_is_initial)) == false
            @test is_clamped(Message(dist2, false, left_is_initial) * Message(dist1, true, right_is_initial)) == false
            @test is_clamped(Message(dist2, true, left_is_initial) * Message(dist1, true, right_is_initial))  == true
        end

        for (left_is_clamped, right_is_clamped) in product(repeated([ true, false ], 2)...)
            @test is_initial(Message(dist1, left_is_clamped, true) * Message(dist2, right_is_clamped, true)) == !(left_is_clamped && right_is_clamped)
            @test is_initial(Message(dist2, left_is_clamped, true) * Message(dist1, right_is_clamped, true)) == !(left_is_clamped && right_is_clamped)
            @test is_initial(Message(dist1, left_is_clamped, false) * Message(dist2, right_is_clamped, false)) == false
            @test is_initial(Message(dist2, left_is_clamped, false) * Message(dist1, right_is_clamped, false)) == false
        end

        @test is_initial(Message(dist1, true, true) * Message(dist2, true, true)) == false
        @test is_initial(Message(dist1, true, true) * Message(dist2, true, false)) == false
        @test is_initial(Message(dist1, true, false) * Message(dist2, true, true)) == false
        @test is_initial(Message(dist1, false, true) * Message(dist2, true, false)) == true
        @test is_initial(Message(dist1, true, false) * Message(dist2, false, true)) == true
        @test is_initial(Message(dist2, true, true) * Message(dist1, true, true)) == false
        @test is_initial(Message(dist2, true, true) * Message(dist1, true, false)) == false
        @test is_initial(Message(dist2, true, false) * Message(dist1, true, true)) == false
        @test is_initial(Message(dist2, false, true) * Message(dist1, true, false)) == true
        @test is_initial(Message(dist2, true, false) * Message(dist1, false, true)) == true
    end

end

end


text:
import Pkg
import UUIDs

@inline function new_compat_entry(current_compat_entry::Vector,
                                  upper_version::VersionNumber)::String
    return new_compat_entry(upper_version)
end

@inline function new_compat_entry(latest_dep_version::VersionNumber)::String
    lower_bound = v"0"
    return new_compat_entry(lower_bound, latest_dep_version)
end

@inline function _extract_lower_bound(current_compat_entry::String)::VersionNumber
    myregex1 = r"([\d]*)\.([\d]*)\.([\d]*)[\s]*-[\s]*"
    myregex2 = r"([\d]*)\.([\d]*)[\s]*-[\s]*"
    myregex3 = r"(\d[\d]*)[\s]*-[\s]*"
    if occursin(myregex1, current_compat_entry)
        m1 = match(myregex1, current_compat_entry)
        return VersionNumber("$(m1[1]).$(m1[2]).$(m1[3])")
    elseif occursin(myregex2, current_compat_entry)
        m2 = match(myregex2, current_compat_entry)
        return VersionNumber("$(m2[1]).$(m2[2])")
    elseif occursin(myregex3, current_compat_entry)
        m3 = match(myregex3, current_compat_entry)
        return VersionNumber("$(m3[1])")
    else
        return v"0"
    end
end

@inline function new_compat_entry(current_compat_entry::String,
                                  upper_version::VersionNumber)::String
    lower_bound = _extract_lower_bound(current_compat_entry)
    return new_compat_entry(lower_bound, upper_version)
end

@inline function _compute_cap_upper_bound(upper_version::VersionNumber)::String
    x = upper_version.major
    y = upper_version.minor
    z = upper_version.patch
    if x == 0
        if y == 0 # x is 0 and y is 0
            return "0.0.$(z)"
        else # x is 0 and y is nonzero
            return "0.$(y)"
        end
    else # x is nonzero
        return "$(x)"
    end
end

@inline function new_compat_entry(lower_bound::VersionNumber,
                                  upper_version::VersionNumber)::String
    a = lower_bound.major
    b = lower_bound.minor
    c = lower_bound.patch
    upper_bound = _compute_cap_upper_bound(upper_version)
    compat = "$(a).$(b).$(c) - $(upper_bound)"
    @assert Pkg.Types.VersionRange(compat) isa Pkg.Types.VersionRange
    return compat
end

@inline function generate_compat_entry(pkg::Package,
                                       dep::Package,
                                       strategy::CapStrategy,
                                       current_compat_entry::Union{Vector, String},
                                       latest_dep_version::VersionNumber,
                                       latest_dep_zero_version::Union{VersionNumber, Nothing})
    if length(current_compat_entry) == 0
        return new_compat_entry(current_compat_entry,
                                latest_dep_version)
    else
        spec = _entry_to_spec(current_compat_entry)
        return generate_compat_entry(pkg,
                                     dep,
                                     strategy,
                                     current_compat_entry,
                                     spec,
                                     latest_dep_version,
                                     latest_dep_zero_version)
    end
end

@inline function _entry_piece_to_ranges(piece::String)::Vector{Pkg.Types.VersionRange}
    try
        return Pkg.Types.semver_spec(piece).ranges
    catch
    end
    r = Pkg.Types.VersionRange(piece)
    return Pkg.Types.VersionRange[r]
end

@inline function _entry_to_spec(current_compat_entry::String)::Pkg.Types.VersionSpec
    try
        return Pkg.Types.semver_spec(current_compat_entry)
    catch
    end
    pieces::Vector{String} = strip.(split(strip(rstrip(lstrip(strip(current_compat_entry), '['), ']')), ','))
    ranges = Vector{Pkg.Types.VersionRange}(undef, 0)
    for piece in pieces
        append!(ranges, _entry_piece_to_ranges(piece))
    end
    spec = Pkg.Types.VersionSpec(ranges)
    return spec
end

@inline function _entry_to_spec(current_compat_entry::Vector{String})::Pkg.Types.VersionSpec
    n = length(current_compat_entry)
    ranges = Vector{Pkg.Types.VersionRange}(undef, n)
    for i = 1:n
        ranges[i] = Pkg.Types.VersionRange(current_compat_entry[i])
    end
    spec = Pkg.Types.VersionSpec(ranges)
end

@inline function generate_compat_entry(pkg::Package,
                                       dep::Package,
                                       strategy::UpperBound,
                                       current_compat_entry::Union{Vector, String},
                                       spec::Pkg.Types.VersionSpec,
                                       latest_dep_version::VersionNumber,
                                       latest_dep_zero_version::Union{VersionNumber, Nothing})
    always_assert(length(current_compat_entry) > 0, "length(current_compat_entry) > 0; pkg=$(pkg), dep=$(dep)")
    if is_unbounded_infinity(spec)
        return new_compat_entry(current_compat_entry, latest_dep_version)
    elseif is_unbounded_bad_zero(spec)
        if latest_dep_zero_version isa Nothing
            return new_compat_entry(current_compat_entry, latest_dep_version)
        else
            return new_compat_entry(current_compat_entry, latest_dep_zero_version)
        end
    else
        # if (next_breaking_release(latest_dep_version) in spec) | (next_breaking_release(latest_dep_zero_version) in spec)
        #     msg = string("This is a bad compat entry because ",
        #                  "it includes future (not-yet-released) ",
        #                  "breaking releases. ",
        #                  "Unfortunately, I am not smart enough to fix this ",
        #                  "compat entry, so I am leaving it unmodified. ",
        #                  "Please manually modify this compat entry.")
        #     @warn(msg,
        #           pkg,
        #           dep,
        #           strategy,
        #           spec,
        #           latest_dep_version,
        #           latest_dep_zero_version,
        #           next_breaking_release(latest_dep_version),
        #           next_breaking_release(latest_dep_zero_version))
        # end
        # if too_many_breaking_zero_versions(spec)
        #     msg = string("This is a bad compat entry because ",
        #                  "it includes an infinite number of breaking ",
        #                  "releases of the form `0.x`.",
        #                  "Unfortunately, I am not smart enough to fix this ",
        #                  "compat entry, so I am leaving it unmodified. ",
        #                  "Please manually modify this compat entry.")
        #     @warn(msg,
        #           pkg,
        #           dep,
        #           strategy,
        #           spec,
        #           latest_dep_version,
        #           latest_dep_zero_version,
        #           next_breaking_release(latest_dep_version),
        #           next_breaking_release(latest_dep_zero_version))
        # end
        return current_compat_entry
    end
end

function too_many_breaking_nonzero_versions(spec::Pkg.Types.VersionSpec, cutoff::Integer = 50)::Bool
    num = 0
    a = typemin(Base.VInt) : Base.VInt(1)       : Base.VInt(200_000)
    b = typemin(Base.VInt) : Base.VInt(100_000) : typemax(Base.VInt)
    for major in Iterators.flatten([a, b])
        if VersionNumber(major, 0, 0) in spec
            num += 1
        end
    end
    if num >= cutoff
        return true
    else
        return false
    end
end

function too_many_breaking_zero_versions(spec::Pkg.Types.VersionSpec, cutoff::Integer = 50)::Bool
    num = 0
    a = typemin(Base.VInt) : Base.VInt(1)       : Base.VInt(200_000)
    b = typemin(Base.VInt) : Base.VInt(100_000) : typemax(Base.VInt)
    for minor in Iterators.flatten([a, b])
        if VersionNumber(0, minor, 0) in spec
            num += 1
        end
    end
    if num >= cutoff
        return true
    else
        return false
    end
end

@inline function includes_infinity(spec::Pkg.Types.VersionSpec)::Bool
    max_component = typemax(Base.VInt)
    x = typemax(VersionNumber)
    y = VersionNumber(max_component, max_component, max_component)
    result = (x in spec) | (y in spec)
    return result
end

@inline function includes_bad_zero(spec::Pkg.Types.VersionSpec)::Bool
    max_component = typemax(Base.VInt)
    a = VersionNumber(0, max_component, 0)
    b = VersionNumber(0, max_component, max_component)
    c = VersionNumber(1, 0, 0)
    has_a_or_b = (a in spec) | (b in spec)
    doesnt_have_c = !(c in spec)
    result = has_a_or_b & doesnt_have_c
    return result
end

@inline function is_unbounded_infinity(spec::Pkg.Types.VersionSpec)::Bool
    # return includes_infinity(spec) | too_many_breaking_nonzero_versions(spec)
    return includes_infinity(spec)
end

@inline function is_unbounded_bad_zero(spec::Pkg.Types.VersionSpec)::Bool
    # return includes_bad_zero(spec) | too_many_breaking_zero_versions(spec)
    return includes_bad_zero(spec)
end

@inline function next_breaking_release(latest_version::VersionNumber)::VersionNumber
    if latest_version < v"1"
        major = 0
        minor = latest_version.minor + 1
    else
        major = latest_version.major + 1
        minor = 0
    end
    return VersionNumber(major, minor, 0)
end

@inline next_breaking_release(::Nothing) = v"0.0.0"

function monotonize!(pkgbounds, compat)
    for (dep, bnd) in compat
        pkgdep = Package(dep)
        if haskey(pkgbounds, pkgdep)
            spec = _entry_to_spec(bnd)
            if !isempty(spec.ranges)
                upper = VersionNumber(spec.ranges[end].upper.t)
                upper == v"0.0.0" && continue
                current_upper = pkgbounds[pkgdep]
                if current_upper === nothing || current_upper == v"0.0.0" || upper < current_upper
                    pkgbounds[pkgdep] = upper
                end
            end
        end
    end
    return pkgbounds
end


text:
# Pressure correlations for PressureDrop package.
# note that all of these work in psia, but the user-facing wrappers for traverses take psig.

#%% Helper functions

"""
`liquidvelocity_superficial(q_o, q_w, id, B_o, B_w)`

Returns superficial liquid velocity, v_sg.

Takes oil rate (q_o, stb/d), water rate (q_w, stb/d), pipe inner diameter (inches), oil volume factor (B_o, dimensionless) and water volume factor (B_w, dimensionless).

Note that this does not account for slip between liquid phases.
"""
function liquidvelocity_superficial(q_o, q_w, id, B_o, B_w)
    A = π * (id/24.0)^2 #convert id in inches to ft

    if q_o > 0
        WOR = q_w / q_o
        return 6.5e-5 * (q_o + q_w) / A * (B_o/(1 + WOR) + B_w * WOR / (1 + WOR))
    else #100% WC
        return 6.5e-5 * q_w * B_w / A
    end
end

"""
`gasvelocity_superficial(q_o, q_w, GLR, R_s, id, B_g)`

Returns superficial liquid velocity, v_sg.

Takes oil rate, (q_o, stb/d), water rate (q_w, stb/d), gas:liquid ratio (scf/stb), solution gas:oil ratio (scf/stb), pipe inner diameter (inches), and gas volume factor (B_g).

Note that this does not account for slip between liquid phases.
"""
function gasvelocity_superficial(q_o, q_w, GLR, R_s, id, B_g)
    A = π * (id/24.0)^2 #convert id in inches to ft

    if q_o > 0
        WOR = q_w / q_o
        return 1.16e-5 * (q_o + q_w) / A * max(GLR - R_s /(1 + WOR), 0) * B_g # max is to prevent edge case where calculated free gas is negative
    else #100% WC
        return 1.16e-5 * q_w * GLR * B_g / A
    end

end

"""
`mixture_properties_simple(q_o, q_w, property_o, property_w)`

Weighted average for mixture properties.

Takes the oil and water rates in equivalent units, as well as their relative properties in equivalent units.

Does not account for oil slip, mixing effects, fluid expansion, behavior of emulsions, etc.
"""
function mixture_properties_simple(q_o, q_w, property_o, property_w)

    return (q_o * property_o + q_w * property_w) / (q_o + q_w)
end

"""
`ChenFrictionFactor(N_Re, id, roughness = 0.01)`

Uses the direct Chen 1979 correlation to determine friction factor, in place of the Colebrook implicit solution.

Takes the dimensionless Reynolds number, pipe inner diameter in **inches**, and roughness in **inches**.

Not intended for Reynolds numbers between 2000-4000.
"""
function ChenFrictionFactor(N_Re, id, roughness = 0.01) #Takacs #returns Economides value * 4

    if N_Re <= 4000 #laminar flow boundary ~2000-2300

        return 64 / N_Re
    else #turbulent flow
        k = roughness/id

        x = -2 * log10(k/3.7065 - 5.0452/N_Re * log10(k^1.1098 / 2.8257 + (7.149/N_Re)^0.8981))

        return 1/x^2
    end
end

"""
`SerghideFrictionFactor(N_Re, id, roughness = 0.01)`

Uses the direct Serghide 1984 correlation to determine friction factor, in place of the Colebrook implicit solution.

Takes the dimensionless Reynolds number, pipe inner diameter in *inches*, and roughness in *inches*.

Not intended for Reynolds numbers between 2000-4000.
"""
function SerghideFrictionFactor(N_Re, id, roughness = 0.01)
    if N_Re <= 4000 #laminar flow boundary ~2000-2300

        return 64 / N_Re
    else #turbulent flow
        k = roughness/id
        A = -2 * log10(k/3.7 + 12/N_Re)
        B = -2 * log10(k/3.7 + 2.51*A/N_Re)
        C = -2 * log10(k/3.7 + 2.51*B/N_Re)

        return (A-((B-A)^2 / (C - (2*B) + A)))^-2
    end
end

#%% Beggs and Brill

"""
`BeggsAndBrillFlowMap(λ_l, N_Fr)`

Beggs and Brill flow pattern as a string ∈ {"segregated", "transition", "distributed", "intermittent"}.

Takes no-slip holdup (λ_l) and mixture Froude number (N_Fr).
"""
function BeggsAndBrillFlowMap(λ_l, N_Fr) #graphical test bypassed in test suite--rerun if modifying this function

    if N_Fr < 316 * λ_l ^ 0.302 && N_Fr < 9.25e-4 * λ_l^-2.468
        return "segregated"
    elseif N_Fr >= 9.25e-4 * λ_l^-2.468 && N_Fr < 0.1 * λ_l^-1.452
        return "transition"
    elseif N_Fr >= 316 * λ_l ^ 0.302 || N_Fr >= 0.5 * λ_l^-6.738
        return "distributed"
    else
        return "intermittent"
    end
end

const BB_coefficients =     (segregated = (a = 0.980, b= 0.4846, c = 0.0868, e = 0.011, f = -3.7680, g = 3.5390, h = -1.6140),
                            intermittent = (a = 0.845, b = 0.5351, c = 0.0173, e = 2.960, f = 0.3050, g = -0.4473, h = 0.0978),
                            distributed = (a = 1.065, b = 0.5824, c = 0.0609),
                            downhill = (e = 4.700, f = -0.3692, g = 0.1244, h = -0.5056) )

"""
`BeggsAndBrillAdjustedLiquidHoldup(flowpattern, λ_l, N_Fr, N_lv, α, inclination, uphill_flow, PayneCorrection = true)`

Helper function for Beggs and Brill. Returns adjusted liquid holdup, ε_l(α), with optional Payne et al correction applied to output.

Takes flow pattern (string ∈ {"segregated", "intermittent", "distributed"}), no-slip holdup (λ_l), Froude number (N_Fr),
liquid velocity number (N_lv), angle from horizontal (α, radians), uphill flow (boolean).
"""
function BeggsAndBrillAdjustedLiquidHoldup(flowpattern, λ_l, N_Fr, N_lv, α, inclination, uphill_flow, PayneCorrection = true)

    if PayneCorrection && uphill_flow
        correctionfactor = 0.924
    elseif PayneCorrection
        correctionfactor = 0.685
    else
        correctionfactor = 1.0
    end

    flow = Symbol(flowpattern)
    a = BB_coefficients[flow][:a]
    b = BB_coefficients[flow][:b]
    c = BB_coefficients[flow][:c]

    ε_l_horizontal = a * λ_l^b / N_Fr^c #liquid holdup assuming horizontal (α = 0 rad)
    ε_l_horizontal = max(ε_l_horizontal, λ_l)

    if α ≈ 0 #horizontal flow
        return ε_l_horizontal
    else #inclined or vertical flow
        if uphill_flow
            if flowpattern == "distributed"
                ψ = 1.0
            else
                e = BB_coefficients[flow][:e]
                f = BB_coefficients[flow][:f]
                g = BB_coefficients[flow][:g]
                h = BB_coefficients[flow][:h]

                C = max( (1 - λ_l) * log(e * λ_l^f * N_lv^g * N_Fr^h), 0)

                if inclination ≈ 0 #vertical flow
                    ψ = 1 + 0.3 * C
                else
                    ψ = 1 + C * (sin(1.8*α) - (1/3) * sin(1.8*α)^3)
                end
            end
        else #downhill flow
            e = BB_coefficients[:downhill][:e]
            f = BB_coefficients[:downhill][:f]
            g = BB_coefficients[:downhill][:g]
            h = BB_coefficients[:downhill][:h]

            C = max( (1 - λ_l) * log(e * λ_l^f * N_lv^g * N_Fr^h), 0)

            if inclination ≈ 0 #vertical flow
                ψ = 1 + 0.3 * C
            else
                ψ = 1 + C * (sin(1.8*α) - (1/3) * sin(1.8*α)^3)
            end
        end

        return ε_l_horizontal * ψ * correctionfactor
    end

end

"""
`BeggsAndBrill(<arguments>)`

Calculates pressure drop for a single pipe segment using Beggs and Brill 1973 method, with optional Payne corrections.

Returns a ΔP in psi.

Doesn't account for oil/water phase slip, but does properly account for inclination.

As of release v0.9, assumes **outlet-defined** models only, i.e. top-down from wellhead; thus, uphill flow corresponds to producers and downhill flow to injectors.

For more information, see *Petroleum Production Systems* by Economides et al., or the the Fekete [reference on pressure drops](http://www.fekete.com/san/webhelp/feketeharmony/harmony_webhelp/content/html_files/reference_material/Calculations_and_Correlations/Pressure_Loss_Calculations.htm).

# Arguments
All arguments take U.S. field units.

- `md`: measured depth of the pipe segment, feet
- `tvd`: true vertical depth, feet
- `inclination`: inclination from vertical, degrees (e.g. vertical => 0)
- `id`: inner diameter of the pipe segment, inches
- `v_sl`: superficial liquid mixture velocity, ft/s
- `v_sg`: superficial gas velocity, ft/s
- `ρ_l`: liquid mixture density, lb/ft³
- `ρ_g`: gas density, lb/ft³
- `σ_l`: liquid/gas interfacial tension, centipoise
- `μ_l`: liquid mixture dynamic viscosity
- `μ_g`: gas dynamic viscosity
- `roughness`: pipe roughness, inches
- `pressure_est`: estimated average pressure of the pipe segment (needed to determine the kinetic effects component of the pressure drop)
- `frictionfactor::Function = SerghideFrictionFactor`: function used to determine the Moody friction factor
- `uphill_flow = true`: indicates uphill or downhill flow. It is assumed that the start of the 1D segment is an outlet and not an inlet
- `PayneCorrection = true`: indicates whether the Payne et al. 1979 corrections should be applied to prevent overprediction of liquid holdup.
"""
function BeggsAndBrill( md, tvd, inclination, id,
                        v_sl, v_sg, ρ_l, ρ_g, σ_l, μ_l, μ_g,
                        roughness, pressure_est, frictionfactor::Function = SerghideFrictionFactor,
                        uphill_flow = true, PayneCorrection = true)

    α = (90 - inclination) * π / 180 #inclination in rad measured from horizontal

    #%% flow pattern and holdup:
    v_m = v_sl + v_sg
    λ_l = v_sl / v_m #no-slip liquid holdup
    N_Fr = 0.373 * v_m^2 / id #mixture Froude number #id is pipe diameter in inches
    N_lv = 1.938 * v_sl * (ρ_l / σ_l)^0.25 #liquid velocity number per Duns & Ros

    flowpattern = BeggsAndBrillFlowMap(λ_l, N_Fr)

    if flowpattern == "transition"
        B = (0.1 * λ_l^-1.4516 - N_Fr) / (0.1 * λ_l^-1.4516 - 9.25e-4 * λ_l^-2.468)
        ε_l_seg = BeggsAndBrillAdjustedLiquidHoldup("segregated", λ_l, N_Fr, N_lv, α, inclination, uphill_flow, PayneCorrection)
        ε_l_int = BeggsAndBrillAdjustedLiquidHoldup("intermittent", λ_l, N_Fr, N_lv, α, inclination, uphill_flow, PayneCorrection)
        ε_l_adj = B * ε_l_seg + (1 - B) * ε_l_int
    else
        ε_l_adj = BeggsAndBrillAdjustedLiquidHoldup(flowpattern, λ_l, N_Fr, N_lv, α, inclination, uphill_flow, PayneCorrection)
    end

    if uphill_flow
        ε_l_adj = max(ε_l_adj, λ_l) #correction to original: for uphill flow, true holdup must by definition be >= no-slip holdup
    end # note that Payne factors reduce the overpredicted liquid holdups from the uncorrected form

    #%% friction factor:
    y = λ_l / ε_l_adj^2
    if 1.0 < y < 1.2
        s = log(2.2y - 1.2) #handle the discontinuity
    else
        ln_y = log(y)
        s = ln_y / (-0.0523 + 3.182 * ln_y - 0.872 * ln_y^2 + 0.01853 * ln_y^4)
    end

    fbyfn = exp(s) #normalizing friction factor f/fₙ

    ρ_ns = ρ_l * λ_l + ρ_g * (1-λ_l) #no-slip density
    μ_ns = μ_l * λ_l + μ_g * (1-λ_l) #no-slip friction in centipoise
    N_Re = 124 * ρ_ns * v_m * id / μ_ns #Reynolds number
    f_n = frictionfactor(N_Re, id, roughness)
    fric = f_n * fbyfn #friction factor

    #%% core calculation:
    ρ_m = ρ_l * ε_l_adj + ρ_g * (1 - ε_l_adj) #mixture density in lb/ft³

    dpdl_el = (1/144.0) * ρ_m
    friction_effect = uphill_flow ? 1 : -1 #note that friction MUST act against the direction of flow
    dpdl_f = friction_effect * 1.294e-3 * fric * (ρ_ns * v_m^2) / id #frictional component
    E_k = 2.16e-4 * fric * (v_m * v_sg * ρ_ns) / pressure_est #kinetic effects; typically negligible

    dp_dl = (dpdl_el * tvd + dpdl_f * md) / (1 - friction_effect*E_k) #assumes friction and kinetic effects both increase pressure in the same 1D direction

    return dp_dl
end #Beggs and Brill

#%% Hagedorn & Brown
"""
`HagedornAndBrownLiquidHoldup(pressure_est, id, v_sl, v_sg, ρ_l, μ_l, σ_l)`

Helper function to determine liquid holdup using the Hagedorn & Brown method.

Does not account for inclination or oil/water slip.
"""
function HagedornAndBrownLiquidHoldup(pressure_est, id, v_sl, v_sg, ρ_l, μ_l, σ_l)

    N_lv = 1.938 * v_sl * (ρ_l / σ_l)^0.25 #liquid velocity number per Duns & Ros
    N_gv = 1.938 * v_sg * (ρ_l / σ_l)^0.25 #gas velocity number per Duns & Ros; yes, use liquid density & viscosity
    N_d = 120.872 * id/12 * (ρ_l / σ_l)^0.5 #pipe diameter number; uses id in ft
    N_l = 0.15726 * μ_l * (1/(ρ_l * σ_l^3))^0.25 #liquid viscosity number

    CN_l = 0.061 * N_l^3 - 0.0929 * N_l^2 + 0.0505 * N_l + 0.0019 #liquid viscosity coefficient * liquid viscosity number

    H = N_lv / N_gv^0.575 * (pressure_est/14.7)^0.1 * CN_l / N_d #holdup correlation group

    ε_l_by_ψ = sqrt((0.0047 + 1123.32 * H + 729489.64 * H^2)/(1 + 1097.1566 * H + 722153.97 * H^2))

    B = N_gv * N_l^0.38 / N_d^2.14
    ψ = (1.0886 - 69.9473*B + 2334.3497*B^2 - 12896.683*B^3)/(1 - 53.4401*B + 1517.9369*B^2 - 8419.8115*B^3) #economides et al 235

    return ψ * ε_l_by_ψ
end

"""
`HagedornAndBrownPressureDrop(pressure_est, id, v_sl, v_sg, ρ_l, ρ_g, μ_l, μ_g, σ_l, id_ft, λ_l, md, tvd, frictionfactor::Function, uphill_flow, roughness)`

Helper function for H&B -- compute H&B pressure drop when bubble flow criteria are not met.
"""
function HagedornAndBrownPressureDrop(pressure_est, id, v_sl, v_sg, ρ_l, ρ_g, μ_l, μ_g, σ_l, id_ft, λ_l, md, tvd, frictionfactor::Function, uphill_flow,  roughness)

    ε_l = HagedornAndBrownLiquidHoldup(pressure_est, id, v_sl, v_sg, ρ_l, μ_l, σ_l)

    if uphill_flow
        ε_l = max(ε_l, λ_l) #correction to original: for uphill flow, true holdup must by definition be >= no-slip holdup
    end

    ρ_m = ρ_l * ε_l + ρ_g * (1 - ε_l) #mixture density in lb/ft³
    massflow = π*(id_ft/2)^2 * (v_sl * ρ_l + v_sg * ρ_g) * 86400 #86400 s/day

    #%% friction factor:
    μ_m = μ_l^ε_l * μ_g^(1-ε_l)
    N_Re = 2.2e-2 * massflow / (id_ft * μ_m)
    fric = frictionfactor(N_Re, id, roughness)/4 #corrected friction factor

    #%% core calculation:
    dpdl_el = 1/144.0 * ρ_m #elevation component
    friction_effect = uphill_flow ? 1 : -1 #note that friction MUST act against the direction of flow
    dpdl_f = friction_effect * 1/144.0 * fric * massflow^2 / (7.413e10 * id_ft^5 *ρ_m) #frictional component: see Economides et al
    #ρ_ns = λ_l * ρ_l + λ_g * ρ_g
    #dpdl_f = 1.294e-3 * fric * ρ_ns^2 * v_m^2 / (ρ_m * id) #Takacs -- takes normal friction factor
    #dpdl_kinetic = 2.16e-4 * ρ_m * v_m * (dvm_dh) #neglected except with high mass flow rates

    dp_dl = dpdl_el * tvd + dpdl_f * md #+ dpdl_kinetic * md

    return dp_dl
end

"""
`GriffithWallisPressureDrop(v_sl, v_sg, v_m, ρ_l, ρ_g, μ_l, id_ft, md, tvd, frictionfactor::Function, uphill_flow, roughness)`

Helper function for H&B correlation -- compute Griffith pressure drop for bubble flow regime.
"""
function GriffithWallisPressureDrop(id, v_sl, v_sg, v_m, ρ_l, ρ_g, μ_l, id_ft, λ_l, md, tvd, frictionfactor::Function, uphill_flow,  roughness)

    v_s = 0.8 #assumed slip velocity of 0.8 ft/s -- probably assumes gas in oil bubbles with no water cut or vice versa?
    ε_l = 1 - 0.5 * (1 + v_m / v_s - sqrt((1 + v_m/v_s)^2 - 4*v_sg/v_s))

    if uphill_flow
        ε_l = max(ε_l, λ_l) #correction to original: for uphill flow, true holdup must by definition be >= no-slip holdup
    end

    ρ_m = ρ_l * ε_l + ρ_g * (1 - ε_l) #mixture density in lb/ft³
    massflow = π*(id_ft/2)^2 * (v_sl * ρ_l + v_sg * ρ_g) * 86400 #86400 s/day

    N_Re = 2.2e-2 * massflow / (id_ft * μ_l) #Reynolds number
    fric = frictionfactor(N_Re, id, roughness)/4 #corrected friction factor

    dpdl_el = 1/144.0 * ρ_m #elevation component
    friction_effect = uphill_flow ? 1 : -1 #note that friction MUST act against the direction of flow
    dpdl_f = friction_effect* 1/144.0 * fric * massflow^2 / (7.413e10 * id_ft^5 * ρ_l * ε_l^2) #frictional component
    dpdl = dpdl_el * tvd + dpdl_f * md

    return dpdl
end

"""
`HagedornAndBrown(<arguments>)`

Calculates pressure drop for a single pipe segment using the Hagedorn & Brown 1965 method (with recent modifications), with optional Griffith and Wallis bubble flow corrections.

Returns a ΔP in psi.

Doesn't account for oil/water phase slip, and does not incorporate flow regimes distinctions outside of in/out of bubble flow. Originally developed for vertical wells.

As of release v0.9, assumes **outlet-defined** models only, i.e. top-down from wellhead; thus, uphill flow corresponds to producers and downhill flow to injectors.

For more information, see *Petroleum Production Systems* by Economides et al., or the the Fekete [reference on pressure drops](http://www.fekete.com/san/webhelp/feketeharmony/harmony_webhelp/content/html_files/reference_material/Calculations_and_Correlations/Pressure_Loss_Calculations.htm).

# Arguments
All arguments take U.S. field units.

- `md`: measured depth of the pipe segment, feet
- `tvd`: true vertical depth, feet
- `inclination`: inclination from vertical, degrees (e.g. vertical => 0)
- `id`: inner diameter of the pipe segment, inches
- `v_sl`: superficial liquid mixture velocity, ft/s
- `v_sg`: superficial gas velocity, ft/s
- `ρ_l`: liquid mixture density, lb/ft³
- `ρ_g`: gas density, lb/ft³
- `σ_l`: liquid/gas interfacial tension, centipoise
- `μ_l`: liquid mixture dynamic viscosity
- `μ_g`: gas dynamic viscosity
- `roughness`: pipe roughness, inches
- `pressure_est`: estimated average pressure of the pipe segment (needed to determine the kinetic effects component of the pressure drop)
- `frictionfactor::Function = SerghideFrictionFactor`: function used to determine the Moody friction factor
- `uphill_flow = true`: indicates uphill or downhill flow. It is assumed that the start of the 1D segment is an outlet and not an inlet
- `GriffithWallisCorrection = true`: indicates whether the Griffith and Wallis 1961 corrections should be applied to prevent overprediction of liquid holdup.
"""
function HagedornAndBrown(md, tvd, inclination, id,
                            v_sl, v_sg, ρ_l, ρ_g, σ_l, μ_l, μ_g,
                            roughness, pressure_est, frictionfactor::Function = SerghideFrictionFactor,
                            uphill_flow = true, GriffithWallisCorrection = true)

    id_ft = id/12

    v_m = v_sl + v_sg

    #%% holdup:
    λ_l = v_sl / v_m
    λ_g = 1 - λ_l

    if GriffithWallisCorrection
        L_B = max(1.071 - 0.2218 * v_m^2 / id, 0.13) #Griffith bubble flow boundary
        if λ_g <  L_B
            dpdl = GriffithWallisPressureDrop(id, v_sl, v_sg, v_m, ρ_l, ρ_g, μ_l, id_ft, λ_l, md, tvd, frictionfactor, uphill_flow, roughness)
        else
            dpdl = HagedornAndBrownPressureDrop(pressure_est, id, v_sl, v_sg, ρ_l, ρ_g, μ_l, μ_g, σ_l, id_ft, λ_l, md, tvd, frictionfactor, uphill_flow, roughness)
        end
    else #no correction
        dpdl = HagedornAndBrownPressureDrop(pressure_est, id, v_sl, v_sg, ρ_l, ρ_g, μ_l, μ_g, σ_l, id_ft, λ_l, md, tvd, frictionfactor, uphill_flow, roughness)
    end

    return dpdl
end #Hagedorn & Brown


text:
"""
    BatchView(data, batchsize; partial=true)
    BatchView(data; batchsize=1, partial=true)

Create a view of the given `data` that represents it as a vector
of batches. Each batch will contain an equal amount of
observations in them. The batch-size
can be specified using the  parameter `batchsize`.
In the case that the size of the dataset is not dividable by the
specified `batchsize`, the remaining observations will
be ignored if `partial=false`. If  `partial=true` instead
the last batch-size can be slightly smaller.

Note that any data access is delayed until `getindex` is called,
and even `getindex` returns an [`ObsView`](@ref)
which in general avoids data movement until [`getobs`](@ref) is
called.

If used as an iterator, the object will iterate over the dataset
once, effectively denoting an epoch.

For `BatchView` to work on some data structure, the type of the
given variable `data` must implement the data container
interface. See [`ObsView`](@ref) for more info.

# Arguments

- **`data`** : The object describing the dataset. Can be of any
    type as long as it implements [`getobs`](@ref) and
    [`numobs`](@ref) (see Details for more information).

- **`batchsize`** : The batch-size of each batch.
    It is the number of observations that each batch must contain
    (except possibly for the last one).

- **`partial`** : If `partial=false` and the number of observations is
    not divisible by the batch-size, then the last mini-batch is dropped.

# Examples

```julia
using MLUtils
X, Y = MLUtils.load_iris()

A = BatchView(X, batchsize=30)
@assert typeof(A) <: BatchView <: AbstractVector
@assert eltype(A) <: SubArray{Float64,2}
@assert length(A) == 5 # Iris has 150 observations
@assert size(A[1]) == (4,30) # Iris has 4 features

# 5 batches of size 30 observations
for x in BatchView(X, batchsize=30)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert numobs(x) === 30
end

# 7 batches of size 20 observations
# Note that the iris dataset has 150 observations,
# which means that with a batchsize of 20, the last
# 10 observations will be ignored
for (x, y) in BatchView((X, Y), batchsize=20, partial=false)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
    @assert numobs(x) == numobs(y) == 20
end

# randomly assign observations to one and only one batch.
for (x, y) in BatchView(shuffleobs((X, Y)), batchsize=20)
    @assert typeof(x) <: SubArray{Float64,2}
    @assert typeof(y) <: SubArray{String,1}
end
```
"""
struct BatchView{TElem,TData} <: AbstractDataContainer
    data::TData
    batchsize::Int
    count::Int
    partial::Bool
    imax::Int
end

function BatchView(data::T; batchsize::Int=1, partial::Bool=true) where {T}
    n = numobs(data)
    if n < batchsize
        @warn "Number of observations less than batch-size, decreasing the batch-size to $n"
        batchsize = n
    end
    E = typeof(obsview(data, 1:batchsize))
    imax = partial ? n : n - batchsize + 1
    count = partial ? ceil(Int, n / batchsize) : floor(Int, n / batchsize)
    BatchView{E,T}(data, batchsize, count, partial, imax)
end

"""
    batchsize(data) -> Int

Return the fixed size of each batch in `data`.
"""
batchsize(A::BatchView) = A.batchsize

numobs(A::BatchView) = A.count
getobs(A::BatchView) = getobs(A.data)
getobs(A::BatchView, i::Int) = getobs(A.data, _batchrange(A, i))

function getobs(A::BatchView, is::AbstractVector)
    obsindices = union((_batchrange(A, i) for i in is)...)::Vector{Int}
    getobs(A.data, obsindices)
end

Base.parent(A::BatchView) = A.data
Base.eltype(::BatchView{Tel}) where Tel = Tel

Base.getindex(A::BatchView, i::Int) = obsview(A.data, _batchrange(A, i))

function Base.getindex(A::BatchView, is::AbstractVector)
    obsindices = union((_batchrange(A, i) for i in is)...)::Vector{Int}
    obsview(A.data, obsindices)
end

obsview(A::BatchView) = A
obsview(A::BatchView, i) = A[i]

# Helper function to translate a batch-index into a range of observations.
function _batchrange(a::BatchView, batchindex::Int)
    (batchindex > a.count || batchindex < 0) && throw(BoundsError())
    startidx = (batchindex - 1) * a.batchsize + 1
    endidx = min(numobs(a.data), startidx + a.batchsize -1)
    return startidx:endidx
end

function Base.showarg(io::IO, A::BatchView, toplevel)
    print(io, "BatchView(")
    Base.showarg(io, parent(A), false)
    print(io, ", ")
    print(io, "batchsize=$(A.batchsize), ")
    print(io, "partial=$(A.partial)")
    print(io, ')')
    toplevel && print(io, " with eltype ", nameof(eltype(A))) # simplify
end

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


text:
using XGPaint
using PyCall
using PyPlot
using WCS
enmap = pyimport("pixell.enmap")

wcs2 = WCSTransform(2;
           cdelt = [1., 1.], #cdelt = [5/60., 5/60.],
           ctype = ["RA---CAR", "DEC--CAR"],
           crpix = [181.0, 91.0],
           crval = [0., 0.])
##
l = collect(1:1000)
cl = l .^ -2.5
shape,wcs = enmap.geometry(shape=(180, 360),
                           res=deg2rad(1),pos=(0,0))
modl = enmap.modlmap(shape, wcs);
cmb = enmap.rand_map(shape, wcs,
   reshape(XGPaint.ellpad(cl),(1,1,length(cl)+1)));
##
imshow(cmb)
gcf()

##

coords = zeros(2,3)
coords = [
   0.0  90.0  180.0;
   0.0  00.0  00.0
]
(world_to_pix(wcs2, coords ))

##
coords = [
   1.0  5.0  180.0;
   1.0  1.0  1.0
]
(pix_to_world(wcs2, coords ))


text:
"""
    eom_chuacircuit_nscroll_sine(u, p, t) -> SVector{3}

Equations of motion for n-scroll chaotic attractors from an adjusted Chua
system [1].

# References
1. Tang, Wallace KS, et al. "Generation of n-scroll attractors via
sine function." IEEE Transactions on Circuits and Systems I:
Fundamental Theory and Applications 48.11 (2001): 1369-1372.

"""
function eom_chuacircuit_nscroll_sine(u, p, t)
    α, β, γ, a, b, c, σx, σy, σz = (p...,)
    x, y, z = (u...,)

    n::Int = c + 1
    if x >= 2*a*c
        fx = (b*pi/2*a)*(x - 2*a*c)
    elseif -2*a*c < x < 2*a*c
        d = ifelse(isodd(n), pi, 0)
        fx = -b*sin((pi*x/2*a) + d)
    elseif x <= -2*a*c
        fx = (b*pi/2*a)*(x + 2*a*c)
    end

    ηx = σx == 0 ? 0 : rand(Normal(0, σx))
    ηy = σy == 0 ? 0 : rand(Normal(0, σy))
    ηz = σz == 0 ? 0 : rand(Normal(0, σz))

    dx = α*(y - fx) + ηx
    dy = x - y + z + ηy
    dz = -β*y - γ*z + ηz
    return SVector{3}(dx, dy, dz)
end

"""
    chua_nscroll_sinefunc(u₀, α, β, γ, a, b, c::Int)

Generate n-scroll chaotic attractors from an adjusted Chua system [1].

# References
1. Tang, Wallace KS, et al. "Generation of n-scroll attractors via
sine function." IEEE Transactions on Circuits and Systems I:
Fundamental Theory and Applications 48.11 (2001): 1369-1372.
"""
function chuacircuit_nscroll_sine(u₀, α, β, γ, a, b, c::Int, σx, σy, σz)
    p = [α, β, γ, a, b, c, σx, σy, σz]
    ContinuousDynamicalSystem(eom_chuacircuit_nscroll_sine, u₀, p)
end
chuacircuit_nscroll_sine(;u₀ = [0.0, 0.0, 0.28695],
        α = 10.814, β = 14, γ = 0, a = 1.3, b = 0.11, c = 2,
        σx = 0.0, σy = 0.0, σz = 0.0) =
    chuacircuit_nscroll_sine(u₀, α, β, γ, a, b, c, σx, σy, σz)


text:
abstract type AbstractTextTransformer <: MMI.Unsupervised end

function MMI.clean!(transformer::AbstractTextTransformer)
    warning = ""
    if transformer.min_doc_freq < 0.0
        warning *= "Need min_doc_freq ≥ 0. Resetting min_doc_freq=0. "
        transformer.min_doc_freq = 0.0
    end

    if transformer.max_doc_freq > 1.0
        warning *= "Need max_doc_freq ≤ 1. Resetting max_doc_freq=1. "
        transformer.max_doc_freq = 1.0
    end

    if transformer.max_doc_freq < transformer.min_doc_freq
        warning *= "max_doc_freq cannot be less than min_doc_freq, resetting to defaults. "
        transformer.min_doc_freq = 0.0
        transformer.max_doc_freq = 1.0
    end
    return warning
end

## General method to fit text transformer models ##
MMI.fit(transformer::AbstractTextTransformer, verbosity::Int, X) =
    _fit(transformer, verbosity, build_corpus(X))

function _fit(transformer::AbstractTextTransformer, verbosity::Int, X::Corpus)
    # process corpus vocab
    update_lexicon!(X)
    dtm_matrix = build_dtm(X)
    n = size(dtm_matrix.dtm, 2) # docs are columns

    # calculate min and max doc freq limits
    if transformer.max_doc_freq < 1 || transformer.min_doc_freq > 0
        high = round(Int, transformer.max_doc_freq * n)
        low = round(Int, transformer.min_doc_freq * n)
        new_doc_term_mat, vocab = limit_features(dtm_matrix, high, low)
    else
        new_doc_term_mat = dtm_matrix.dtm
        vocab = dtm_matrix.terms
    end

    # calculate IDF
    idf = compute_idf(transformer.smooth_idf, new_doc_term_mat)

    # prepare result
    fitresult = get_result(transformer, idf, vocab, new_doc_term_mat)
    cache = nothing

    return fitresult, cache, NamedTuple()
end

## General method to transform using text transformer models ##
MMI.transform(transformer::AbstractTextTransformer, result, v) =
    _transform(transformer, result, build_corpus(v))


text:
"""
```
compute_moments(m::AbstractModel, percent::Float64 = 0.90; verbose::Symbol=:none)
```

Computes prior and posterior parameter moments. Tabulates prior mean, posterior mean, and
bands in various LaTeX tables stored `tablespath(m)`.

### Arguments
- `m`: the model object
- `percent`: the percentage of the mass of draws from Metropolis-Hastings included between
  the bands displayed in output tables.
"""
function compute_moments(m::AbstractModel, percent::AbstractFloat = 0.90; verbose::Symbol=:none)

    # Read in the matrix of parameter draws from metropolis-hastings
    filename = rawpath(m, "estimate", "mhsave.h5")

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        @printf "Reading parameter draws from %s\n" filename
    end

    param_draws = []
    try
        h5open(filename, "r") do f
            param_draws = read(f, "mhparams")
        end
    catch
        @printf "Could not open file %s\n" filename
        return
    end

    n_draws = size(param_draws,1)

    # Save mean parameter vector
    save_mean_parameters(m, param_draws)

    # Produce TeX table of moments
    make_moment_tables(m, param_draws, percent, verbose=verbose)
end

"""
```
save_mean_parameters{T<:AbstractFloat}(m::AbstractModel, draws::Matrix{T})
```

Computes and saves the posterior mean of the parameters.

### Arguments
- `m`
- `draws`: n_draws x n_parameters matrix holding the posterior draws from
  Metropolis-Hastings
"""
function save_mean_parameters{T<:AbstractFloat}(m::AbstractModel, draws::Matrix{T})
    post_means = mean(draws,1)'
    filename = workpath(m, "estimate", "paramsmean.h5")
    h5open(filename, "w") do f
        f["params"] = post_means
    end
end

"""
```
make_moment_tables{T<:AbstractFloat}(m::AbstractModel, draws::Matrix{T},
    percent::AbstractFloat; verbose::Symbol = :none)
```

Tabulates parameter moments in 3 LaTeX tables:

1. For MAIN parameters, a list of prior means, prior standard deviations, posterior means,
   and 90% bands for posterior draws

2. For LESS IMPORTANT parameters, a list of the prior means, prior standard deviations,
   posterior means and 90% bands for posterior draws.

3. A list of prior means and posterior means

### Arguments
- `draws`: [n_draws x n_parameters] matrix holding the posterior draws from
  Metropolis-Hastings from save/mhsave.h5
- `percent`: the mass of observations we want; 0 <= percent <= 1
"""
function make_moment_tables{T<:AbstractFloat}(m::AbstractModel, draws::Matrix{T},
                                              percent::AbstractFloat;
                                              verbose::Symbol = :none)

    # STEP 1: Extract moments of prior distribution from m.parameters
    n_params    = length(m.parameters)
    prior_means = zeros(n_params,1)
    prior_std   = zeros(n_params,1)

    # Track prior distributions
    prior_dist  = Vector{AbstractString}(n_params)
    distid(::Distributions.Beta)    = "B"
    distid(::Distributions.Gamma)   = "G"
    distid(::Distributions.Normal)  = "N"
    distid(::DSGE.RootInverseGamma) = "IG"

    for (i,k) in enumerate(m.keys)

        if i > n_params
            continue
        end

        param = getindex(m,i)

        # If param is fixed, we simply set mean to its value and std to 0
        if param.fixed
            prior_means[i]  = param.value
            prior_std[i] = 0
            prior_dist[i] = "-"
        else
            pri = get(param.prior)
            α, β = moments(pri)
            prior_means[i] = α
            prior_std[i] = β
            prior_dist[i] = distid(pri)
        end
    end

    # STEP 2: Compute moments and `percent' bands from parameter draws

    # Posterior mean for each
    post_means = mean(draws,1)'

    # Covariance: Note that this has already been computed and saved
    # in $workpath(m)/parameter_covariance.h5.
    # parameter θ_sig = cov(θ, mean=θ_hat')

    # Bands for each
    post_bands = []
    try
        post_bands = find_density_bands(draws, percent, minimize=true)
    catch
        println("percent must be between 0 and 1")
        return -1
    end

    # We need the transpose
    post_bands = post_bands'

    # STEP 3: Create variables for writing to output table

    # Two-row column names. First row is multicolumn, where entries (i,str) are `i` columns
    # with content `str`.
    colnames0 = [(1,""), (3,"Prior"),(3,"Posterior")]
    colnames = ["Parameter", "Type", "Mean", "SD", "Mean",
                "$(100*percent)\\% {\\tiny Lower Band}",
                "$(100*percent)\\% {\\tiny Upper Band}"]

    # prior type, mean, std dev; posterior mean, 90% lower band, 90% upper bands
    outmat = [prior_dist prior_means prior_std post_means post_bands]

    # prior mean and posterior mean (n_params x 2)
    outmat2 = [prior_means post_means]

    # STEP 4: WRITE TABLES

    # 4a. Write to Table 1: prior mean, std dev, posterior mean and bands for IMPORTANT
    #     parameters

    # Open and start the TeX file
    moments0_out = tablespath(m, "estimate", "moments0.tex")
    moments0_fid = open(moments0_out, "w")

    write_table_preamble(moments0_fid)

    @printf moments0_fid "\\caption{Parameter Estimates}\n"
    @printf moments0_fid "\\vspace*{.5cm}\n"
    @printf moments0_fid "{\\small \n"
    @printf moments0_fid "\\begin{tabular}{lcccccc}\n"
    @printf moments0_fid "\\hline\n"

    # Column names
    @printf moments0_fid "\\multicolumn{%d}{c}{%s}" colnames0[1][1] colnames0[1][2]
    for col in colnames0[2:end]
        @printf moments0_fid " & \\multicolumn{%d}{c}{%s}" col[1] col[2]
    end
    @printf moments0_fid "\\\\ \n"
    @printf moments0_fid "%s" colnames[1]
    for col in colnames[2:end]
        @printf moments0_fid " & %s" col
    end
    @printf moments0_fid "\\\\ \n"
    @printf moments0_fid "\\cmidrule(lr){1-1}\\cmidrule(lr){2-4}\\cmidrule(lr){5-7}\n"

    # Keep track of indices for important parameters
    important_para = []

    for (index, param) in enumerate(m.parameters)

        if (!ismatch(r"rho_", param.tex_label) &&
            !ismatch(r"zeta_", param.tex_label) &&
            !ismatch(r"psi_", param.tex_label) &&
            !ismatch(r"nu_l", param.tex_label) &&
            !ismatch(r"pi\^\*", param.tex_label) &&
            !ismatch(r"sigma_{pi}\^\*",param.tex_label) &&
            (!ismatch(r"pistar", param.tex_label)))
            continue
        end

        if ismatch(r"rho_chi",param.tex_label)
            continue
        end

        #Print the parameter name and values in outmat
        @printf moments0_fid "\$\%4.99s\$ & " param.tex_label
        @printf moments0_fid "%s & %8.3f & %8.3f & %8.3f & %8.3f & %8.3f \\\\\n" outmat[index,:]...

        important_para = [important_para; index]
    end

    # Close the file
    write_table_postamble(moments0_fid;
        note="{\\footnotesize Note: For Inverse Gamma (IG) prior mean and SD, \$\\tau^2\$ and \$\\nu\$
        reported.}",
        small=true)

    # 4b. Write to Table 2: Prior mean, std dev and posterior mean, bands for other params
    moments_table_out = tablespath(m, "estimate", "moments1.tex")
    moments_table_fid = open(moments_table_out, "w")

    write_table_preamble(moments_table_fid)

    @printf moments_table_fid "\\caption{Parameter Estimates}\n"
    @printf moments_table_fid "\\vspace*{.2cm}\n"
    @printf moments_table_fid "{\\small \n"
    @printf moments_table_fid "\\begin{tabular}{lcccccc}\n"
    @printf moments_table_fid "\\hline\n"

    # Column names
    for col in colnames0
        @printf moments_table_fid "\\multicolumn{%d}{c}{%s} &" col[1] col[2]
    end
    @printf moments_table_fid "\\\n"
    for col in colnames
        @printf moments_table_fid "%s & " col
    end
    @printf moments_table_fid "\\\n"
    @printf moments_table_fid "\\cmidrule(lr){1-1}\\cmidrule(lr){2-4}\\cmidrule(lr){5-7}\n"

    # Counter for parameters to track length of table and number of tables in excess of
    # default (1)
    other_para  = 1
    table_count = 0

    for (index, param) in enumerate(m.parameters)

        if in(index, important_para)
            continue
        end

        # Make a new table if the current one is too large
        if other_para % 25 == 0 && index ≠ length(m.parameters)

            # Close the file
            write_table_postamble(moments_table_fid;
                note="{\\footnotesize Note: For Inverse Gamma (IG) prior mean and SD, \$\\tau^2\$ and \$\\nu\$
                reported.}",
                small=true)

            # Update table counter
            table_count += 1

            # Start the new file
            filename = @sprintf "moments%d.tex" table_count
            moments_table_out = tablespath(m,"estimate",filename)
            moments_table_fid = open(moments_table_out,"w")

            write_table_preamble(moments_table_fid)
            @printf moments_table_fid "\\caption{Parameter Estimates}\n"
            @printf moments_table_fid "\\vspace*{.2cm}\n"
            @printf moments_table_fid "{\\small \n"
            @printf moments_table_fid "\\begin{tabular}{lcccccc}\n"
            @printf moments_table_fid "\\hline\n"

            # Column names
            @printf moments_table_fid "\\multicolumn{%d}{c}{%s}" colnames0[1][1] colnames0[1][2]
            for col in colnames0[2:end]
                @printf moments_table_fid " & \\multicolumn{%d}{c}{%s}" col[1] col[2]
            end
            @printf moments_table_fid "\\\\ \n"
            @printf moments_table_fid "%s" colnames[1]
            for col in colnames[2:end]
                @printf moments_table_fid " & %s" col
            end
            @printf moments_table_fid "\\\\ \n"
            @printf moments_table_fid "\\cmidrule(lr){1-1}\\cmidrule(lr){2-4}\\cmidrule(lr){5-7}\n"
        end

        #Print the parameter name and values in outmat
        @printf moments_table_fid "\$\%4.99s\$ & " param.tex_label
        @printf moments_table_fid "%s & %8.3f & %8.3f & %8.3f & %8.3f & %8.3f \\\\\n" outmat[index,:]...

        other_para += 1
    end

    write_table_postamble(moments_table_fid; small=true)

    # 4c. Write to Table 5: Prior mean and posterior mean for all parameters

    # Keep track of how many tables we've made
    table_count = 0

    # Open the TeX file and set up the heading
    prioPostMean_out = tablespath(m,"estimate", "moments_prioPostMean.tex")
    prioPostMean_fid = open(prioPostMean_out,"w")

    write_table_preamble(prioPostMean_fid)
    @printf prioPostMean_fid "\\caption{Parameter Estimates: Prior and Posterior Mean}\n"
    @printf prioPostMean_fid "\\vspace*{.5cm}\n"
    @printf prioPostMean_fid "\\begin{tabular}{ccc}\\hline \n"
    @printf prioPostMean_fid " Parameter & Prior & Posterior  \\tabularnewline \\hline\n"

    # Write out the results
    for (index, param) in enumerate(m.parameters)

        if index % 40 == 0 && index ≠ length(m.parameters)

            # Close the old file
            write_table_postamble(prioPostMean_fid)

            # Generate the new filename
            table_count += 1
            filename = @sprintf "moments_prioPostMean_%d.tex" table_count

            # Open a new file and start the next table
            prioPostMean_out = tablespath(m,"estimate",filename)
            prioPostMean_fid = open(prioPostMean_out,"w")
            write_table_preamble(prioPostMean_fid)
            @printf prioPostMean_fid "\\caption{Parameter Estimates: Prior and Posterior Mean}\n"
            @printf prioPostMean_fid "\\vspace*{.5cm}\n"
            @printf prioPostMean_fid "\\begin{tabular}{ccc}\\hline \n"
            @printf prioPostMean_fid " Parameter & Prior & Posterior  \\tabularnewline \\hline\n"
        end

        @printf prioPostMean_fid "\$\%4.99s\$ & " param.tex_label

        val = outmat2[index,:]
        @printf prioPostMean_fid "\%8.3f & \%8.3f \\\\\n" val[1] val[2]
    end

    write_table_postamble(prioPostMean_fid)

    if VERBOSITY[verbose] >= VERBOSITY[:low]
        @printf "Tables are saved as %s.\n" tablespath(m, "estimate", "*.tex")
    end
end

function write_table_preamble(fid::IOStream)
    @printf fid "\\documentclass[12pt]{article}\n"
    @printf fid "\\usepackage{booktabs}\n"
    @printf fid "\\begin{document}\n"
    @printf fid "\\pagestyle{empty}\n"
    @printf fid "\\begin{table}[h]\n"
    @printf fid "\\centering\n"
end

# `small`: Whether to print an additional curly bracket after "\end{tabular}" (necessary if
# the table is enclosed by "\small{}")

function write_table_postamble(fid::IOStream; note::AbstractString="", small::Bool=false)
    @printf fid "\\\\ \\\hline\n"

    if small
        @printf fid "\\end{tabular}}\n"
    else
        @printf fid "\\end{tabular}\n"
    end

    if !isempty(note)
        @printf fid "%s\n" note
    end

    @printf fid "\\end{table}\n"
    @printf fid "\\end{document}"
    close(fid)
end

"""
```
find_density_bands(draws::Matrix, percent::AbstractFloat; minimize::Bool=true)
```

Returns a `2` x `cols(draws)` matrix `bands` such that `percent` of the mass of `draws[:,i]`
is above `bands[1,i]` and below `bands[2,i]`.

### Arguments
- `draws`: Matrix of parameter draws (from Metropolis-Hastings, for example)
- `percent`: percent of data within bands (e.g. .9 to get 90% of mass within bands)

### Optional Arguments
- `minimize`: if `true`, choose shortest interval, otherwise just chop off lowest and
  highest (percent/2)
"""
function find_density_bands(draws::Matrix, percent::AbstractFloat; minimize::Bool=true)

    if !(0 < percent < 1)
        error("percent must be between 0 and 1")
    end

    n_draws, n_draw_dimensions = size(draws)
    band = zeros(2, n_draw_dimensions)
    n_in_band  = round(Int, percent * n_draws)

    for i in 1:n_draw_dimensions

        # Sort response for parameter i such that 1st element is largest
        draw_variable_i = draws[:,i]
        sort!(draw_variable_i, rev=true)

        # Search to find the interval containing the minimum # of observations
        # comprising `percent` of the mass
        if minimize

            upper_index = 1
            done        = 0
            j           = 2
            minwidth = draw_variable_i[1] - draw_variable_i[n_in_band]

            while j <= (n_draws - n_in_band + 1)

                newwidth = draw_variable_i[j] - draw_variable_i[j + n_in_band - 1]

                if newwidth < minwidth
                    upper_index = j
                    minwidth = newwidth
                end

                j += 1
            end

        else
            upper_index = n_draws - nwidth - floor(.5*n_draws-n_in_band)
        end

        band[2,i] = draw_variable_i[upper_index]
        band[1,i] = draw_variable_i[upper_index + n_in_band - 1]
    end

    return band
end


text:
using Documenter
using LiterateTest
using Literate

rm(joinpath(@__DIR__, "src", "examples"), force = true, recursive = true)
Literate.markdown(
    joinpath(@__DIR__, "..", "examples", "tests.jl"),
    joinpath(@__DIR__, "src", "examples");
    documenter = true,
)

makedocs(
    sitename = "LiterateTest",
    format = Documenter.HTML(),
    modules = [LiterateTest],
    pages = [
        "Home" => "index.md",
        "Testing LiterateTest.jl" => "examples/tests.md",
    ],
)

deploydocs(; repo = "github.com/tkf/LiterateTest.jl", push_preview = true)


text:
type OnePointCrossover{R} <: Crossover
  stage::AbstractString
  rep::R # Is this really necessary?
  rate::Float

  OnePointCrossover(stage::AbstractString, rep::R, rate::Float) =
    new(stage, rep, rate)
end

"""
Used for composing one point crossover operators.
"""
type OnePointCrossoverDefinition <: CrossoverDefinition
  stage::AbstractString
  rate::Float

  OnePointCrossoverDefinition() = new("", 0.7)
  OnePointCrossoverDefinition(r::Float) = new("", r)
  OnePointCrossoverDefinition(s::AbstractString) = new(s, 0.7)
  OnePointCrossoverDefinition(s::AbstractString, r::Float) = new(s, r)
end

"""
Composes a one-point crossover operator from a provided definition, for a given
species.
"""
function compose!(c::OnePointCrossoverDefinition, sp::Species)
  c.stage = c.stage == "" ? genotype(sp).label : c.stage
  r = sp.stages[c.stage].representation
  OnePointCrossover{typeof(r)}(c.stage, r, c.rate)
end

"""
TODO: Description of one-point crossover.
"""
one_point() = one_point(0.7)
one_point(rate::Float) = OnePointCrossoverDefinition(rate)
one_point(stage::AbstractString) = OnePointCrossoverDefinition(stage, 0.7)
one_point(stage::AbstractString, rate::Float) =
  OnePointCrossoverDefinition(stage, rate)
function one_point(def::Function)
  d = one_point()
  def(d)
  d
end

num_inputs(o::OnePointCrossover) = 2
num_outputs(o::OnePointCrossover) = 2

function operate!{T}(o::OnePointCrossover, outputs::Vector{IndividualStage{Vector{T}}}, x::Vector{T}, y::Vector{T})
  # Enforce the crossover rate.
  rand() > o.rate && return

  # Cache the length of x and y.
  lx = length(x)
  ly = length(y)

  # Ensure that the inputs are greater than length 1 (should we do this dynamically?).
  (lx <= 1 || ly <= 1)

  # Calculate the crossover point, split A and B into four substrings
  # and combine those substrings to form two children.
  c = rand(2:(min(lx, ly) - 1))

  # Generate the output chromosomes.
  t = splice!(x, (c + 1):ly, y[(c+1):end])
  splice!(y, (c + 1):lx, t)

  return
end


text:
###Target functions
function sincos!(r::Matrix,paras::Vector,t::AbstractVector)
    p1 = 2pi*paras[1]
    p2 = 2pi*paras[2]
    @simd for i=1:length(t)
        @inbounds r[i,1] = sin(p1*t[i])
    end
    @simd for i=1:length(t)
        @inbounds r[i,2] = cos(p2*t[i])
    end
    r
end

sincos(paras::Vector,t::AbstractVector) = sincos!(zeros(eltype(t),length(t),2),paras,t)

@inline function _commonsincoscomponents(t::AbstractVector,measurementparas::Vector,variance::Vector,paraminit...)
    modelparameters = parameters([:a,:b],paraminit...)
    noisemodel = noise(:gaussian,variance)
    measurementdata = data(:array,t,applynoise!(noisemodel,sincos(measurementparas,t)))
    modelparameters,measurementdata,noisemodel
end

###Create a target function model
function sincosmodel(t::AbstractVector,measurementparas::Vector,variance::Vector,paraminit...)
    model(:target,_commonsincoscomponents(t,measurementparas,variance,paraminit...)...,sincos,t;name="Sine-Cosine Model")
end

###Create a target function model
function sincosmodel!(t::AbstractVector,measurementparas::Vector,variance::Vector,paraminit...)
    model(:target!,_commonsincoscomponents(t,measurementparas,variance,paraminit...)...,sincos!,t;name="Sine-Cosine In-Place Model")
end


text:
# Testing sparseUtils
using jInv.Utils
using Test

# test sdiag
a = randn(13)
@test a == diag(sdiag(a))


text:
#!/usr/bin/env julia
# del2z <delta.z@aliyun.com>

module GloVe

using Flux
using Flux: glorot_normal, ADAM, OneHotMatrix, OneHotVector, onehot, onehotbatch, onecold, logitbinarycrossentropy, throttle
using JSON
using DataStructures
using StatsBase

const config = JSON.parsefile(joinpath(@__DIR__, "word2vec.json"))

include("prepro.jl")

function dataset end
function model end
function loss end

end


text:
# zadatak 8.2
include("labyrinth.jl")
include("zadatak82BFS.jl")

# kreiramo graf
nbRows = 8;
nbColumns = 8;
G = emptyLabyrinth(nbRows, nbColumns);

# kreiramo lavirint
defineNodesType!(G, 1:nbRows*nbColumns, 'P');
defineNodesType!(G, [5 9 10 11 13 15 21 23 24 26 28 29 31 34 44 45 46 47 50 55 58 60 61 63], 'Z');

# definisemo izlaze
Izlazi=[1 8 62 64]

# trazimo put od 57 do najblizeg izlaza
path = BFS!(G, 57, Izlazi)
println("Path = $path")

printLabyrinthPath(G, path)


text:
function test_parse(input, args... ; broken=false)
    arg = (:expr, args...)
    if broken
        @test_broken texparse(input) == manual_texexpr(arg)
    else
        @test texparse(input) == manual_texexpr(arg)
    end
end

@testset "Parser" begin
    @testset "Accent" begin
        # First char is the combining accent
        test_parse(
            raw"\dot{a}",
            (:combining_accent, '̇', 'a')
        )
        test_parse(
            raw"\vec{x}",
            (:combining_accent, '⃗', 'x')
        )
    end

    @testset "Delimiter" begin
        test_parse(
            raw"\left( x \right)",
            (:delimited, '(', 'x', ')')
        )

        test_parse(
            raw"\left[ y \right)",
            (:delimited, '[', 'y', ')')
        )

        test_parse(
            raw"\left( a + b \right)",
            (:delimited,
                '(',
                (:group, 'a', (:spaced, '+'), 'b'),
                ')'
            )
        )

        @test_throws TeXParseError texparse(raw"\left( x")
        @test_throws TeXParseError texparse(raw"x \right)")

        # TODO Recursive delim
    end

    @testset "Fonts" begin
        test_parse(raw"\mathrm{u}", (:font, :rm, (:char, 'u')))
    end

    @testset "Fraction" begin
        test_parse(raw"\frac{1}{n}", (:frac, (:digit, '1'), 'n'))
    end

    @testset "Function" begin
        # Check that braces are not added to the function name
        test_parse(raw"\sin{x}", (:function, "sin"), 'x')
        test_parse(
            raw"\exp(x)",
            (:function, "exp"), '(', 'x', ')'
        )
        test_parse(
            raw"\tan\alpha",
            (:function, "tan"), (:symbol, 'α')
        )
    end

    @testset "Group" begin
        test_parse(raw"{x}", 'x')
        test_parse(raw"{fgh}", (:group, 'f', 'g', 'h'))
        test_parse(
            raw"{a{b{cd}}}",
            (:group, 'a', (:group, 'b', (:group, 'c', 'd')))
        )

        @test_throws TeXParseError texparse(raw"{v")
        @test_throws TeXParseError texparse(raw"w}")
    end

    @testset "Integral" begin
        test_parse(
            raw"\int",
            (:integral, (:symbol, '∫'), nothing, nothing)
        )
        test_parse(
            raw"\int_a^b",
            (:integral, (:symbol, '∫'), 'a', 'b')
        )
    end

    @testset "Overunder" begin
        test_parse(
            raw"\sum",
            (:underover, (:symbol, '∑'), nothing, nothing)
        )
        test_parse(
            raw"\sum_{k=0}^n",
            (:underover,
                (:symbol, '∑'),
                (:group, 'k', (:spaced, (:symbol, '=')), (:digit, '0')),
                'n'
            )
        )
        test_parse(
            raw"\lim_x",
            (:underover,
                (:function, "lim"),
                'x',
                nothing
            )
        )
    end

    @testset "Square root" begin
        test_parse(raw"\sqrt{x}", (:sqrt, 'x'))
        test_parse(
            raw"\sqrt{abc}",
            (:sqrt, (:group, 'a', 'b', 'c'))
        )
    end

    @testset "Space" begin
        test_parse(raw"\quad", (:space, 1))
        test_parse(raw"\qquad", (:space, 2))
        test_parse(raw"~", (:space, 0.33333))
    end

    @testset "Spaced symbol" begin
        test_parse(raw"=", (:spaced, (:symbol, '=')))
        test_parse(
            raw"\rightarrow",
            (:spaced, (:symbol, '→'))
        )
        # Hyphen must be replaced by a minus sign
        test_parse(raw"-", (:spaced, (:symbol, '−')))
    end

    @testset "Subscript and superscript" begin
        @test texparse(raw"a^2_3") == texparse("a_3^2")
    end

    @testset "Symbol" begin
        for (char, sym) in zip(split("ϕ φ Φ"), split(raw"\phi \varphi \Phi"))
            test_parse(sym, (:symbol, first(char)))
            @test texparse(char) == texparse(sym)
        end

        # Check interaction with decoration
        test_parse(
            raw"ω_k",
            (:decorated, (:symbol, 'ω'), 'k', nothing)
        )
    end
end