text:
@doc raw"""
    HyperbolicDiffusionEquations3D

The linear hyperbolic diffusion equations in three space dimensions.
A description of this system can be found in Sec. 2.5 of the book "I Do Like CFD, Too: Vol 1".
The book is freely available at http://www.cfdbooks.com/ and further analysis can be found in
the paper by Nishikawa [DOI: 10.1016/j.jcp.2007.07.029](https://doi.org/10.1016/j.jcp.2007.07.029)
"""
struct HyperbolicDiffusionEquations3D <: AbstractHyperbolicDiffusionEquations{3, 4}
  Lr::Float64
  Tr::Float64
  nu::Float64
  resid_tol::Float64
end

function HyperbolicDiffusionEquations3D()
  # diffusion coefficient
  nu = parameter("nu", 1.0)
  # relaxation length scale
  Lr = parameter("Lr", 1.0/(2.0*pi))
  # relaxation time
  Tr = Lr*Lr/nu
  # stopping tolerance for the pseudotime "steady-state"
  resid_tol = parameter("resid_tol", 1e-12)
  HyperbolicDiffusionEquations3D(Lr, Tr, nu, resid_tol)
end

get_name(::HyperbolicDiffusionEquations3D) = "HyperbolicDiffusionEquations3D"
varnames_cons(::HyperbolicDiffusionEquations3D) = @SVector ["phi", "q1", "q2", "q3"]
varnames_prim(::HyperbolicDiffusionEquations3D) = @SVector ["phi", "q1", "q2", "q3"]
default_analysis_quantities(::HyperbolicDiffusionEquations3D) = (:l2_error, :linf_error, :residual)

# Set initial conditions at physical location `x` for pseudo-time `t`
function initial_conditions_poisson_periodic(x, t, equation::HyperbolicDiffusionEquations3D)
  # elliptic equation: -νΔϕ = f
  # depending on initial constant state, c, for phi this converges to the solution ϕ + c
  if iszero(t)
    phi = 0.0
    q1  = 0.0
    q2  = 0.0
    q3  = 0.0
  else
    phi =          sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * x[3])
    q1  = 2 * pi * cos(2 * pi * x[1]) * sin(2 * pi * x[2]) * sin(2 * pi * x[3])
    q2  = 2 * pi * sin(2 * pi * x[1]) * cos(2 * pi * x[2]) * sin(2 * pi * x[3])
    q3  = 2 * pi * sin(2 * pi * x[1]) * sin(2 * pi * x[2]) * cos(2 * pi * x[3])
  end
  return @SVector [phi, q1, q2, q3]
end

function initial_conditions_poisson_nonperiodic(x, t, equation::HyperbolicDiffusionEquations3D)
  # elliptic equation: -νΔϕ = f
  if t == 0.0
    phi = 1.0
    q1  = 1.0
    q2  = 1.0
    q3  = 1.0
  else
    phi =  2.0 *      cos(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) + 2.0 # ϕ
    q1  = -2.0 * pi * sin(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3])   # ϕ_x
    q2  =  4.0 * pi * cos(pi * x[1]) * cos(2.0 * pi * x[2]) * sin(2.0 * pi * x[3])   # ϕ_y
    q3  =  4.0 * pi * cos(pi * x[1]) * sin(2.0 * pi * x[2]) * cos(2.0 * pi * x[3])   # ϕ_z
  end
  return @SVector [phi, q1, q2, q3]
end

function initial_conditions_eoc_test_coupled_euler_gravity(x, t, equation::HyperbolicDiffusionEquations3D)

  # Determine phi_x, phi_y
  G = 1.0 # gravitational constant
  C_grav = -4 * G / (3 * pi) # "3" is the number of spatial dimensions  # 2D: -2.0*G/pi
  A = 0.1 # perturbation coefficient must match Euler setup
  rho1 = A * sin(pi * (x[1] + x[2] + x[3] - t))
  # intialize with ansatz of gravity potential
  phi = C_grav * rho1
  q1  = C_grav * A * pi * cos(pi*(x[1] + x[2] + x[3] - t)) # = gravity acceleration in x-direction
  q2  = q1                                                 # = gravity acceleration in y-direction
  q3  = q1                                                 # = gravity acceleration in z-direction

  return @SVector [phi, q1, q2, q3]
end

function initial_conditions_sedov_self_gravity(x, t, equation::HyperbolicDiffusionEquations3D)
  # for now just use constant initial condition for sedov blast wave (can likely be improved)

  phi = 0.0
  q1  = 0.0
  q2  = 0.0
  q3  = 0.0
  return @SVector [phi, q1, q2, q3]
end

# Apply boundary conditions
function boundary_conditions_sedov_self_gravity(u_inner, orientation, direction, x, t,
                                                surface_flux_function,
                                                equation::HyperbolicDiffusionEquations3D)
  u_boundary = initial_conditions_sedov_self_gravity(x, t, equation)

  # Calculate boundary flux
  if direction in (2, 4, 6) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

function boundary_conditions_poisson_nonperiodic(u_inner, orientation, direction, x, t,
                                                 surface_flux_function,
                                                 equation::HyperbolicDiffusionEquations3D)
  # elliptic equation: -νΔϕ = f
  phi =  2.0 *      cos(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3]) + 2.0 # ϕ
  q1  = -2.0 * pi * sin(pi * x[1]) * sin(2.0 * pi * x[2]) * sin(2.0 * pi * x[3])   # ϕ_x
  q2  =  4.0 * pi * cos(pi * x[1]) * cos(2.0 * pi * x[2]) * sin(2.0 * pi * x[3])   # ϕ_y
  q3  =  4.0 * pi * cos(pi * x[1]) * sin(2.0 * pi * x[2]) * cos(2.0 * pi * x[3])   # ϕ_z
  u_boundary = @SVector [phi, q1, q2, q3]

  # Calculate boundary flux
  if direction in (2, 4, 6) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equation)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equation)
  end

  return flux
end

# Apply source terms
function source_terms_poisson_periodic(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations3D)
  # elliptic equation: -νΔϕ = f
  # analytical solution: phi = sin(2πx)*sin(2πy)*sin(2πz) and f = -12νπ^2 sin(2πx)*sin(2πy)*sin(2πz)
  inv_Tr = inv(equation.Tr)
  C = -12.0*equation.nu*pi*pi

  for k in 1:n_nodes, j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, k, element_id]
    x2 = x[2, i, j, k, element_id]
    x3 = x[3, i, j, k, element_id]
    tmp1 = sin(2 * pi * x1)
    tmp2 = sin(2 * pi * x2)
    tmp3 = sin(2 * pi * x3)
    ut[1, i, j, k, element_id] -= C*tmp1*tmp2*tmp3
    ut[2, i, j, k, element_id] -= inv_Tr * u[2, i, j, k, element_id]
    ut[3, i, j, k, element_id] -= inv_Tr * u[3, i, j, k, element_id]
    ut[4, i, j, k, element_id] -= inv_Tr * u[4, i, j, k, element_id]
  end

  return nothing
end

function source_terms_poisson_nonperiodic(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations3D)
  # elliptic equation: -νΔϕ = f
  # analytical solution: ϕ = 2 cos(πx)sin(2πy)sin(2πz) + 2 and f = 18 π^2 cos(πx)sin(2πy)sin(2πz)
  inv_Tr = inv(equation.Tr)

  for k in 1:n_nodes, j in 1:n_nodes, i in 1:n_nodes
    x1 = x[1, i, j, k, element_id]
    x2 = x[2, i, j, k, element_id]
    x3 = x[3, i, j, k, element_id]
    ut[1, i, j, k, element_id] += 18 * pi^2 * cos(pi*x1) * sin(2.0*pi*x2) * sin(2.0*pi*x3)
    ut[2, i, j, k, element_id] -= inv_Tr * u[2, i, j, k, element_id]
    ut[3, i, j, k, element_id] -= inv_Tr * u[3, i, j, k, element_id]
    ut[4, i, j, k, element_id] -= inv_Tr * u[4, i, j, k, element_id]
  end

  return nothing
end

function source_terms_harmonic(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations3D)
  # harmonic solution ϕ = (sinh(πx)sin(πy) + sinh(πy)sin(πx))/sinh(π), so f = 0
  inv_Tr = inv(equation.Tr)

  for k in 1:n_nodes, j in 1:n_nodes, i in 1:n_nodes
    ut[2, i, j, k, element_id] -= inv_Tr * u[2, i, j, k, element_id]
    ut[3, i, j, k, element_id] -= inv_Tr * u[3, i, j, k, element_id]
    ut[4, i, j, k, element_id] -= inv_Tr * u[4, i, j, k, element_id]
  end

  return nothing
end

# The coupled EOC test does not require additional sources
function source_terms_eoc_test_coupled_euler_gravity(ut, u, x, element_id, t, n_nodes, equation::HyperbolicDiffusionEquations3D)
  return source_terms_harmonic(ut, u, x, element_id, t, n_nodes, equation)
end

# Calculate 1D flux in for a single point
@inline function calcflux(u, orientation, equation::HyperbolicDiffusionEquations3D)
  phi, q1, q2, q3 = u

  if orientation == 1
    f1 = -equation.nu*q1
    f2 = -phi/equation.Tr
    f3 = zero(phi)
    f4 = zero(phi)
  elseif orientation == 2
    f1 = -equation.nu*q2
    f2 = zero(phi)
    f3 = -phi/equation.Tr
    f4 = zero(phi)
  else
    f1 = -equation.nu*q3
    f2 = zero(phi)
    f3 = zero(phi)
    f4 = -phi/equation.Tr
  end

  return SVector(f1, f2, f3, f4)
end

@inline function flux_lax_friedrichs(u_ll, u_rr, orientation, equation::HyperbolicDiffusionEquations3D)
  # Obtain left and right fluxes
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  λ_max = sqrt(equation.nu / equation.Tr)

  return 0.5 * (f_ll + f_rr - λ_max * (u_rr - u_ll))
end

@inline function flux_upwind(u_ll, u_rr, orientation, equation::HyperbolicDiffusionEquations3D)
  # Obtain left and right fluxes
  phi_ll, q1_ll, q2_ll, q3_ll = u_ll
  phi_rr, q1_rr, q2_rr, q3_rr = u_rr
  f_ll = calcflux(u_ll, orientation, equation)
  f_rr = calcflux(u_rr, orientation, equation)

  # this is an optimized version of the application of the upwind dissipation matrix:
  #   dissipation = 0.5*R_n*|Λ|*inv(R_n)[[u]]
  λ_max = sqrt(equation.nu/equation.Tr)
  f1 = 1/2 * (f_ll[1] + f_rr[1]) - 1/2 * λ_max * (phi_rr - phi_ll)
  if orientation == 1 # x-direction
    f2 = 1/2 * (f_ll[2] + f_rr[2]) - 1/2 * λ_max * (q1_rr - q1_ll)
    f3 = 1/2 * (f_ll[3] + f_rr[3])
    f4 = 1/2 * (f_ll[4] + f_rr[4])
  elseif orientation == 2 # y-direction
    f2 = 1/2 * (f_ll[2] + f_rr[2])
    f3 = 1/2 * (f_ll[3] + f_rr[3]) - 1/2 * λ_max * (q2_rr - q2_ll)
    f4 = 1/2 * (f_ll[4] + f_rr[4])
  else # y-direction
    f2 = 1/2 * (f_ll[2] + f_rr[2])
    f3 = 1/2 * (f_ll[3] + f_rr[3])
    f4 = 1/2 * (f_ll[4] + f_rr[4]) - 1/2 * λ_max * (q3_rr - q3_ll)
  end

  return SVector(f1, f2, f3, f4)
end

# Determine maximum stable time step based on polynomial degree and CFL number
function calc_max_dt(u, element_id, invjacobian, cfl,
                     equation::HyperbolicDiffusionEquations3D, dg)
  λ_max = sqrt(equation.nu / equation.Tr)
  dt = cfl * 2 / (nnodes(dg) * invjacobian * λ_max)

  return dt
end

# Convert conservative variables to primitive
@inline cons2prim(u, equation::HyperbolicDiffusionEquations3D) = u

# Convert conservative variables to entropy found in I Do Like CFD, Too, Vol. 1
@inline function cons2entropy(u, equation::HyperbolicDiffusionEquations3D)
  phi, q1, q2, q3 = u
  w1 = phi
  w2 = equation.Lr^2 * q1
  w3 = equation.Lr^2 * q2
  w4 = equation.Lr^2 * q3

  return SVector(w1, w2, w3, w4)
end

# Calculate entropy for a conservative state `u` (here: same as total energy)
@inline entropy(u, equation::HyperbolicDiffusionEquations3D) = energy_total(u, equation)

# Calculate total energy for a conservative state `u`
@inline function energy_total(u, equation::HyperbolicDiffusionEquations3D)
  # energy function as found in equation (2.5.12) in the book "I Do Like CFD, Vol. 1"
  phi, q1, q2, q3 = u
  return 0.5 * (phi^2 + equation.Lr^2 * (q1^2 + q2^2 + q3^2))
end


text:
using AutoDiff
using Scalar
using FactCheck

facts( "Univariate" ) do

    a = Dual( 1.0, 2.0 )
    @fact typeof( a ) <: Dual --> true
    @fact isa( a, Dual{Float64,Float64}) --> true
    @fact re( a ) --> 1.0
    @fact du( a ) --> 2.0

    b = Dual( 3.0, 4.0 )
    @fact typeof( a ) <: Dual --> true
    @fact isa( b, Dual{Float64,Float64}) --> true
    @fact re( b ) --> 3.0
    @fact du( b ) --> 4.0

    test_binary_op( op::Symbol, re_expected::Float64, du_expected::Float64 ) = begin
        context( string( op ) ) do
            c = eval(op)( a, b )
            @fact typeof( c ) <: Dual --> true
            @fact isa( c, Dual{Float64, Float64}) --> true
            @fact re( c ) --> re_expected
            @fact du( c ) --> du_expected
        end
    end

    test_binary_op( :+, re( a )+re( b ), du( a )+du( b ))
    test_binary_op( :-, re( a )-re( b ), du( a )-du( b ))
    test_binary_op( :*, re( a )*re( b ), re( a )*du( b )+re( b )*du( a ))
    test_binary_op( :/, re( a )/re( b ), du( a )/re( b )-du( b )*re( a )/re( b )^2 )
    test_binary_op( :^, re( a )^re( b ), du( a )*re( b )*re( a )^( re( b )-1 )+du( b )*re( a )^re( b )*log( re( a )))

    test_unary_op( op::Symbol, re_expected::Float64, du_expected::Float64 ) = begin
        context( string( op )) do
            c = eval( op )( a )
            @fact typeof( c ) <: Dual --> true
            @fact isa( c, Dual{Float64,Float64}) --> true
            @fact re( c ) --> re_expected
            @fact du( c ) --> du_expected
        end
    end

    test_unary_op( :cos, cos( re( a )), -sin( re( a ))*du( a ))
    test_unary_op( :sin, sin( re( a )), cos( re( a ))*du( a ))
    test_unary_op( :tan, tan( re( a )), (tan( re( a ))^2+1)*du( a ))

end

facts( "Multivariate" ) do
end

facts( "Sparse Multivariate" ) do
end

ext: jl

max_stars_repo_name: sabjohnso/AutoDiff

max_issues_repo_name: sabjohnso/AutoDiff

max_forks_repo_name: sabjohnso/AutoDiff


text:
module TestAsType

using Test
using ArbFloats

af0 = zero(ArbFloat)
af1 = one(ArbFloat)
af2 = af1 + af1
afhalf = af1 / 2
afNaN = ArbFloat{precision(ArbFloat)}(NaN)
afInf = ArbFloat(Inf)

@test isnan(afNaN)
@test isinf(afInf)
@test isfinite(af1)
@test !isfinite(afInf)

@test af1 === af1
@test af1 == af1
@test af1 != af2

@test typeof(ArbFloat{53}(af2)) == ArbFloat{53}
@test typeof(af2) == ArbFloat{precision(ArbFloat)}

@test string(ArbFloat{24}(pi)) == "3.141593"

@test promote(ArbFloat(4), 8//2) == (ArbFloat(4), ArbFloat(4))

end # module


text:
using BenchmarkTools

#=
println("preparing Eirene")
include(joinpath(@__DIR__, "eirene-comparison.jl"))
e = run(EireneComparison.suite())
println("running Eirene")
BenchmarkTools.save("results-eirene.json", e)

println("preparing Matrix")
include(joinpath(@__DIR__, "matrix-reduction.jl"))
m = run(MatrixReduction.suite())
println("running Matrix")
BenchmarkTools.save("results-matrix.json", m)

println("preparing Scaling")
include(joinpath(@__DIR__, "scaling.jl"))
run_gcycle_and_random16_scaling()
run_subset_scaling()
=#


text:
module ClickDetection

# using AxisArrays
using IntervalSets
using FileIO
using LibSndFile
using SampledSignals
using Statistics
using DataFrames
using Dates

include("detector.jl")
include("io.jl")

export read_as_axisarray,
    AbstractClick,
    ClickPointer,
    Click,
    left,
    right,
    series,
    interval,
    samplerate,
    samples,
    times,
    teager,
    merge_intervals,
    detect_clicks,
    peak_lead_ratio,
    as_dataframe

end # module


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

"""
    ArrayShape{T,N} <: AbstractValueShape

Describes the shape of `N`-dimensional arrays of type `T` and a given size.

Constructor:

    ArrayShape{T}(dims::NTuple{N,Integer}) where {T,N}
    ArrayShape{T}(dims::Integer...) where {T}

e.g.

    shape = ArrayShape{Real}(2, 3)

See also the documentation of [`AbstractValueShape`](@ref).
"""
struct ArrayShape{T,N} <: AbstractValueShape
    dims::NTuple{N,Int}
end

export ArrayShape

ArrayShape{T}(dims::NTuple{N,Integer}) where {T,N} = ArrayShape{T,N}(map(Int, dims))
ArrayShape{T}(dims::Integer...) where {T} = ArrayShape{T}(dims)

@inline Base.size(shape::ArrayShape) = shape.dims
Base.length(shape::ArrayShape) = prod(size(shape))

import Base.<=
@inline <=(a::ArrayShape{T,N}, b::ArrayShape{U,N}) where {T,U,N} = T<:U && size(a) == size(b)

@inline _valshapeoftype(T::Type{<:AbstractArray}) = throw(ArgumentError("Type $T does not have a fixed shape"))

@inline default_unshaped_eltype(shape::ArrayShape{T}) where {T} =
    default_unshaped_eltype(_valshapeoftype(T))

@inline shaped_type(shape::ArrayShape{T,N}, ::Type{U}) where {T,N,U<:Real} =
    Array{shaped_type(_valshapeoftype(T),U),N}

@inline function valshape(x::AbstractArray{T}) where T
    _valshapeoftype(T) # ensure T has a fixed shape
    ArrayShape{T}(size(x))
end

@inline function valshape(x::AbstractArray{T,0}) where T
    _valshapeoftype(T) # ensure T has a fixed shape
    ScalarShape{T}()
end

# Possible extension: valshape(x::AbstractArrayOfSimilarArrays{...})

totalndof(shape::ArrayShape{T}) where{T} =
    prod(size(shape)) * totalndof(_valshapeoftype(T))

(shape::ArrayShape{T,N})(::UndefInitializer) where {T,N} = Array{default_datatype(T),N}(undef, size(shape)...)

function _check_unshaped_compat(A::AbstractArray{T,N}, shape::ArrayShape{U,N}) where {T<:Real,U<:Real,N}
    Telem = eltype(A)
    Telem <: U || throw(ArgumentError("Element type $Telem of array not compatible with element type $U of given shape"))
    size(A) == size(shape) || throw(ArgumentError("Size of array differs from size of given shape"))
end

function unshaped(A::AbstractArray{T,1}, shape::ArrayShape{U,1}) where {T<:Real,U<:Real}
    _check_unshaped_compat(A, shape)
    A
end

function unshaped(A::AbstractArray{T,N}, shape::ArrayShape{U,N}) where {T<:Real,U<:Real,N}
    _check_unshaped_compat(A, shape)
    reshape(view(A, ntuple(_ -> :, Val(N))...), prod(size(A)))
end

function unshaped(A::Base.ReshapedArray{T,N,<:AbstractArray{T,1}}, shape::ArrayShape{U,N}) where {T<:Real,U<:Real,N}
    _check_unshaped_compat(A, shape)
    parent(A)
end

replace_const_shapes(f::Function, shape::ArrayShape) = shape

@inline function _apply_shape_to_data(shape::ArrayShape{<:Real,1}, data::AbstractVector{<:Real})
    @boundscheck _checkcompat(shape, data)
    data
end

const ArrayAccessor{T,N} = ValueAccessor{ArrayShape{T,N}} where {T,N}

const RealScalarOrVectorAccessor = ValueAccessor{<:Union{ScalarShape{<:Real},ArrayShape{<:Real,1}}}

Base.@propagate_inbounds vs_getindex(data::AbstractVector{<:Real}, va::ArrayAccessor{<:Real}) = copy(view(data, va))

@static if VERSION < v"1.4"
    # To avoid ambiguity with Julia v1.0 (and v1.1 to v1.3?)
    Base.@propagate_inbounds vs_getindex(data::AbstractVector{<:Real}, va::ArrayAccessor{<:Real,1}) = copy(view(data, va))
end

Base.@propagate_inbounds function vs_getindex(
    data::AbstractArray{<:Real,N},
    idxs::Vararg{RealScalarOrVectorAccessor,N}
) where N
    idxs_mapped = map(view_idxs, axes(data), idxs)
    getindex(data, idxs_mapped...)
end

Base.@propagate_inbounds vs_unsafe_view(data::AbstractVector{<:Real}, va::ArrayAccessor{<:Real,1}) =
    Base.unsafe_view(data, view_idxs(axes(data, 1), va))

Base.@propagate_inbounds vs_unsafe_view(data::AbstractVector{<:Real}, va::ArrayAccessor{<:Real,N}) where {N} =
    reshape(Base.unsafe_view(data, view_idxs(axes(data, 1), va)), size(va.shape)...)

Base.@propagate_inbounds function vs_unsafe_view(
    data::AbstractArray{<:Real,N},
    idxs::Vararg{RealScalarOrVectorAccessor,N}
) where N
    idxs_mapped = map(view_idxs, axes(data), idxs)
    Base.view(data, idxs_mapped...)
end

Base.@propagate_inbounds vs_setindex!(data::AbstractVector{<:Real}, v, va::ArrayAccessor{<:Real}) =
    setindex!(data, v, view_idxs(axes(data, 1), va))

@static if VERSION < v"1.4"
    # To avoid ambiguity with Julia v1.0 (and v1.1 to v1.3?)
    Base.@propagate_inbounds vs_setindex!(data::AbstractVector{<:Real}, v, va::ArrayAccessor{<:Real,1}) =
        setindex!(data, v, view_idxs(axes(data, 1), va))
end

Base.@propagate_inbounds function vs_setindex!(
    data::AbstractArray{<:Real,N},
    v,
    idxs::Vararg{RealScalarOrVectorAccessor,N}
) where N
    idxs_mapped = map(view_idxs, axes(data), idxs)
    setindex!(data, v, idxs_mapped...)
end

Base.@propagate_inbounds function _bcasted_view(data::AbstractVectorOfSimilarVectors{<:Real}, va::ArrayAccessor{T,1}) where {T}
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    fpview = view(flat_data, idxs, :)
    VectorOfSimilarVectors(fpview)
end

Base.@propagate_inbounds function _bcasted_view(data::AbstractVectorOfSimilarVectors{<:Real}, va::ArrayAccessor{T,N}) where {T,N}
    flat_data = flatview(data)
    idxs = view_idxs(axes(flat_data, 1), va)
    fpview = view(flat_data, idxs, :)
    VectorOfSimilarArrays(reshape(fpview, size(va.shape)..., :))
end

Base.copy(instance::VSBroadcasted2{1,typeof(getindex),AbstractVectorOfSimilarVectors{<:Real},Ref{<:ArrayAccessor}}) =
    copy(_bcasted_view(instance.args[1], instance.args[2][]))

Base.copy(instance::VSBroadcasted2{1,typeof(view),AbstractVectorOfSimilarVectors{<:Real},Ref{<:ArrayAccessor}}) =
    _bcasted_view(instance.args[1], instance.args[2][])

function _bcasted_view_unchanged(data::AbstractArray{<:AbstractVector{T}}, shape::ArrayShape{U,1}) where {T<:Real,U>:T}
    _checkcompat_inner(shape, data)
    data
end

# Specialize (::ArrayShape{T,1}).(::AbstractArray{<:AbstractVector{<:Real}}):
Base.copy(instance::VSBroadcasted1{N,ArrayShape{T,1},AbstractArray{<:AbstractVector{<:Real},N}}) where {T,N} =
    _bcasted_view_unchanged(instance.args[1], instance.f)

@inline _bcasted_unshaped(A::AbstractArrayOfSimilarArrays{<:Real,1}) = A

# Specialize unshaped.(::AbstractArray{<:AbstractVector{<:Real}}):
Base.copy(instance::VSBroadcasted1{N,typeof(unshaped),AbstractArray{<:AbstractVector{<:Real},N}}) where N =
    _bcasted_unshaped(instance.args[1])

# TODO: Add support for StaticArray.

# Possible extension: variable/flexible array shapes?


text:
using DataFrames, Gadfly

set_default_plot_size(9inch, 3inch)

shapes = [Shape.square, Shape.diamond, Shape.cross, Shape.xcross, Shape.utriangle, Shape.dtriangle,
    Shape.star1, Shape.star2, Shape.hexagon, Shape.octagon, Shape.hline, Shape.vline,
    Shape.ltriangle, Shape.rtriangle]

D = DataFrame(x=[1:5; 1:4; 1:3; 1:2], y=[5:-1:1; 4:-1:1; 3:-1:1; 2:-1:1], g=string.(1:14))

function plotf(coord::Coord.cartesian)
    plot(D, coord,
        layer(x=:x, y=:y, Theme(default_color="black", point_size=2px)),
        layer(x=:x, y=:y, shape=:g, color=:g, Geom.point),
    Theme(point_shapes=shapes, key_position=:none, point_size=4mm,
     plot_padding=[0mm], background_color="white") )
end

coords = [Coord.cartesian(xflip=xf, yflip=yf, xmin=0, xmax=6, ymin=0, ymax=6)
    for (xf, yf) in zip([false, true, false], [false, false, true]) ]
plots = plotf.(coords)
hstack(plots...)


text:
# module DoubleDouble

# export Double, Single, double
import Base.convert, Base.*, Base.+, Base.-, Base./, Base.sqrt, Base.<, Base.rem,
Base.abs, Base.rand, Base.promote_rule, Base.one, Base.zero, Base.ones,
Base.zeros, Base.show

const BitsFloat=Union{BigFloat,Float16,Float32,Float64} # Floating point BitTypes AbstractFloat

abstract type AbstractDouble <: Number end

# a Single is a wrapper for an ordinary floating point type such that arithmetic operations will return Doubles
struct Single{T} <: AbstractDouble
    hi::T
end

# In a Double, hi uses the full mantissa, and abs(lo) <= 0.5eps(hi)
struct Double{T} <: AbstractDouble
    hi::T
    lo::T
end

Double(x::T) where {T<:BitsFloat} = Double(x,zero(T))  # {T<:BitsFloat}

const half64 = 1.34217729e8
const half32 = 4097f0
const half16 = Float16(33.0)
const halfBig = 3.402823669209384634633746074317682114570000000000000000000000000000000000000000e+38
# 6.805647338418769269267492148635364229120000000000000000000000000000000000000000e38

# round floats to half-precision
# TODO: fix overflow for large values
halfprec(x::Float64) = (p = x*half64; (x-p)+p) # signif(x,26,2) for 26 is 6.7108865e7, this seems like 27
halfprec(x::Float32) = (p = x*half32; (x-p)+p) # float32(signif(x,12,2))
halfprec(x::Float16) = (p = x*half16; (x-p)+p) # float16(signif(x,5,2))
halfprec(x::BigFloat) = (p = x*halfBig; (x-p)+p) # BigFloat(signif(x,128,2))

function splitprec(x::BitsFloat)
    h = halfprec(x)
    h, x-h
end

# Dodano
one(x::Double) = oftype(x,1.0)
one(::Type{T}) where {T<:Double} = convert(T,1.0)
zero(x::Double) = oftype(x,0.0)
zero(::Type{T}) where {T<:Double} = convert(T,0.0)
show(x::Double) = show(STDOUT::IO, [x.hi,x.lo])
ones(T::Double, dims...) = fill!(Array(T, dims...), (one)(T))
zeros(T::Double, dims...) = fill!(Array(T, dims...), (zero)(T))

convert(::Type{Double{T}}, x::Int64) where T = Double(Float64(x))
double(x::Int64) = convert(Double{Float64},x)
Double(x::Int64) = double(x)
promote_rule(::Type{Double{T}},::Type{Int64}) where {T<:BitsFloat} =Double{T}
promote_rule(::Type{Int64},::Type{Double{T}}) where {T<:BitsFloat}=Double{T}

## conversion and promotion
convert(::Type{Single{T}}, x::T) where {T<:BitsFloat} = Single(x)
convert(::Type{Double{T}}, x::T) where {T<:BitsFloat} = Double(x)

convert(::Type{Double{T}}, x::Single{T}) where {T<:BitsFloat} = Double(x.hi)
convert(::Type{Single{T}}, x::Double{T}) where {T<:BitsFloat} = Single(x.hi)

# convert(::Type{T}, x::AbstractDouble{T}) where {T<:BitsFloat} = x.hi

convert(::Type{Single{T}}, x::Single{T}) where {T<:BitsFloat} = x # needed because Double <: FloatingPoint
convert(::Type{Double{T}}, x::Double{T}) where {T<:BitsFloat} = x # needed because Double <: FloatingPoint

convert(::Type{Single{T}}, x::AbstractFloat) where {T<:BitsFloat} = Single(convert(T,x))

function convert(::Type{Double{T}}, x::AbstractFloat) where {T<:BitsFloat}
    z = convert(T,x)
    Double(z,convert(T,x-z))
end

convert(::Type{BigFloat}, x::Single{T}) where {T<:BitsFloat} = big(x.hi)
convert(::Type{BigFloat}, x::Double{T}) where {T<:BitsFloat} = big(x.hi) + big(x.lo)
convert(::Type{Float64}, x::Double{T}) where {T<:BitsFloat} = Float64(x.hi) + Float64(x.lo)

promote_rule(::Type{Single{T}}, ::Type{T}) where {T<:BitsFloat} = Single{T}
promote_rule(::Type{Double{T}}, ::Type{T}) where {T<:BitsFloat} = Double{T}
promote_rule(::Type{Double{T}}, ::Type{Single{T}}) where {T<:BitsFloat} = Double{T}

# promote_rule{T<:BitsFloat}(::Type{AbstractDouble{T}}, ::Type{BigFloat}) = BigFloat  !!
promote_rule(::Type{Irrational{s}}, ::Type{Single{T}}) where {s,T<:BitsFloat} = Double{BigFloat}

double(x::BitsFloat) = Double(x)
# "Normalise" doubles to ensure abs(lo) <= 0.5eps(hi)
# assumes abs(u) > abs(v): if not, use Single + Single
# could be moved to the constructor?
function double(u::T,v::T) where {T<:BitsFloat}
    w = u + v
    Double(w,(u-w) + v)
end
double(x::BigFloat) = convert(Double{BigFloat},x)
double(x::Irrational{S}) where S = convert(Double{Float64},x)

# <

function <(x::Double{T},y::Double{T}) where T
    x.hi+x.lo < y.hi+y.lo ? true : false
end

# add12
function +(x::Single{T},y::Single{T}) where T
    abs(x.hi) > abs(y.hi) ? double(x.hi,y.hi) : double(y.hi,x.hi)
end

# Dekker add2
function +(x::Double{T}, y::Double{T}) where T
    r = x.hi + y.hi
    s = abs(x.hi) > abs(y.hi) ? (((x.hi - r) + y.hi) + y.lo) + x.lo : (((y.hi - r) + x.hi) + x.lo) + y.lo
    double(r,s)
end

# add122
function +(x::Single{T}, y::Double{T}) where T
    r = x.hi + y.hi
    s = abs(x.hi) > abs(y.hi) ? ((x.hi - r) + y.hi) + y.lo : ((y.hi - r) + x.hi) + y.lo
    double(r,s)
end
+(x::Double{T}, y::Single{T}) where T = y + x

-(x::Double{T}) where {T<:BitsFloat} = Double(-x.hi,-x.lo)

function -(x::Double{T}, y::Double{T}) where T
    r = x.hi - y.hi
    s = abs(x.hi) > abs(y.hi) ? (((x.hi - r) - y.hi) - y.lo) + x.lo : (((-y.hi - r) + x.hi) + x.lo) - y.lo
    double(r,s)
end

# Dekker mul12
function *(x::Single{T},y::Single{T}) where T
    hx,lx = splitprec(x.hi)
    hy,ly = splitprec(y.hi)
    z = x.hi*y.hi
    Double(z, ((hx*hy-z) + hx*ly + lx*hy) + lx*ly)
end

# Dekker mul2
function *(x::Double{T}, y::Double{T}) where T
    c = Single(x.hi) * Single(y.hi)
    cc = (x.hi * y.lo + x.lo* y.hi) + c.lo
    double(c.hi, cc)
end

# Dekker div2
function /(x::Double{T}, y::Double{T}) where T
    c = x.hi / y.hi
    u = Single(c) * Single(y.hi)
    cc = ((((x.hi - u.hi) - u.lo) + x.lo) - c*y.lo)/y.hi
    double(c,cc)
end

# Dekker sqrt2
function sqrt(x::Double{T}) where T
    if x.hi <= 0
        throw(DomainError("sqrt will only return a complex result if called with a complex argument."))
    end
    c = sqrt(x.hi)
    u = Single(c)*Single(c)
    cc = (((x.hi - u.hi) - u.lo) + x.lo)*map(typeof(x.hi),0.5)/c
    double(c,cc)
end

rem(x::Double{T},d::Real) where T = double(rem(x.hi,d),rem(x.lo,d))
abs(x::Double{T}) where T = x.hi>0 ? x : -x

# random numbers using full Uint64 range (respectively, UInt32, UInt16 and UInt128)
function rand(::Type{Double{Float64}})
    u = rand(UInt64)
    f = Float64(u)
    uf = UInt64(f)
    ur = uf > u ? uf-u : u-uf
    Double(5.421010862427522e-20*f, 5.421010862427522e-20*Float64(ur))
end

function rand(::Type{Double{Float32}})
    u = rand(UInt32)
    f = Float32(u)
    uf = UInt32(f)
    ur = uf > u ? uf-u : u-uf
    Double(2.3283064f-10*f, 2.3283064f-10*Float32(ur))
end

function rand(::Type{Double{Float16}})
    u = rand(UInt16)
    f = Float16(u)
    uf = UInt16(f)
    ur = uf > u ? uf-u : u-uf
    Double(Float16(1.526e-5)*f, Float16(1.526e-5)*Float16(ur))
end

function rand(::Type{Double{BigFloat}})
    u = rand(UInt128)
    f = BigFloat(u)
    uf = UInt128(f)
    ur = uf > u ? uf-u : u-uf
    Double(2.938735877055718769921841343055614194546663891930218803771879265696043148636818e-39*f, 2.938735877055718769921841343055614194546663891930218803771879265696043148636818e-39*BigFloat(ur))
end

# calculate constants from big numbers
macro twofloat_const_frombig(sym)
    esym = esc(sym)
    qsym = esc(Expr(:quote, sym))
    bigval = @eval big($sym)
    quote
        Base.convert(::Type{Double{Float64}}, ::Irrational{$qsym}) = $(convert(Double{Float64}, bigval))
        Base.convert(::Type{Double{Float32}}, ::Irrational{$qsym}) = $(convert(Double{Float32}, bigval))
        Base.convert(::Type{Double{Float16}}, ::Irrational{$qsym}) = $(convert(Double{Float16}, bigval))
        Base.convert(::Type{Double{BigFloat}}, ::Irrational{$qsym}) = $(convert(Double{BigFloat}, bigval))
    end
end

#=
@twofloat_const_frombig π
@twofloat_const_frombig MathConstants.e
@twofloat_const_frombig MathConstants.γ
@twofloat_const_frombig MathConstants.catalan
@twofloat_const_frombig φ
=#

# end #module


text:
export ir2jl
function ir2jl(directory::String, o::String)
    decls = IR.Decl[]
    for each in readdir(directory)
        each = joinpath(directory, each)
        isfile(each) || continue
        endswith(each, ".mlfso") || continue

        d = open(each) do f
            JSON.Parser.parse(read(f, String))
        end
        solvedIRs = fromVec(Vector{IR.Decl}, d).value
        append!(decls, solvedIRs)
    end
    ir = IR.ELet(decls, IR.Expr(LineNumberNode(1),nothing, IR.EInt(0, 8)))
    jl_exp = irToJulia(ir, LineNumberNode(1))
    open(o, "w") do ff
        println(ff, Expr(:let, Expr(:block), jl_exp))
    end
end


text:
using DataFrames
using LinearAlgebra
using TraitSimulation
using Random

#this VarianceComponent type stores A, B , CholA and CholB so we don't have to compute the cholesky decomposition inside the loop

struct VarianceComponent
	A::Matrix{Float64} # n_traits by n_traits
	B::Matrix{Float64} # n_people by n_people
	CholA::Cholesky{Float64,Array{Float64,2}} # cholesky decomposition of A
	CholB::Cholesky{Float64,Array{Float64,2}} # cholesky decomposition of B

	function VarianceComponent(A, B) #inner constructor given A, B
		return(new(A, B, cholesky(A), cholesky(B))) # stores these values (this is helpful so we don't have it inside the loop)
	end
end

#single LMM trait with given evaluated matrix of Variance components
function LMM_trait_simulation(mu, vc::Matrix{T}) where T
	n_people = size(mu)[1]
	n_traits = size(mu)[2]
	#preallocate memory for the returned dataframe simulated_trait
	simulated_trait = zeros(n_people, n_traits)
	z = Matrix{Float64}(undef, n_people, n_traits)

	#for a single evaluated matrix as the specified covariance matrix instead of a Variancewe do not need to call the
	chol_Σ = cholesky(vc)
	#generate from standard normal
	randn!(z)

	# we want to solve u then v to get the first variane component, v.
	#first matrix vector multiplication using cholesky decomposition

	#need to find which will be CholA, CholB
	lmul!(chol_Σ.L, z)

	simulated_trait += z

	#for each trait
	simulated_trait += mu

	out = DataFrame(simulated_trait)

	out = names!(out, [Symbol("trait$i") for i in 1:n_traits])

	return out
end

function append_terms!(AB, summand)
	# elements in args are symbols,
	A_esc = esc(summand.args[2])
	B_esc = esc(summand.args[3])
	push!(AB.args, :(VarianceComponent($A_esc, $B_esc)))
end

"""
this vc macro allows us to create a vector of VarianceComponent objects for simulation
"""

macro vc(expression)
	n = length(expression.args)
	# AB is an empty vector of variance components list of symbols
	AB = :(VarianceComponent[])
	if expression.args[1] != :+ #if first argument is not plus (only one vc)
		summand = expression
		append_terms!(AB, summand)
	else #MULTIPLE VARIANCE COMPONENTS if the first argument is a plus (Sigma is a sum multiple variance components)
		for i in 2:n
			summand = expression.args[i]
			append_terms!(AB, summand)
		end
	end
	return(:($AB)) # change this to return a vector of VarianceComponent objects
end

"""
this is a test for vcobjtuple that is compatible with VarianceComponentModels.jl
"""
function  vcobjtuple(vcobject::Vector{VarianceComponent})
	m = length(vcobject)
	d = size(vcobject[1].A, 1)
	n = size(vcobject[1].B, 1)

	Σ = ntuple(x -> zeros(d, d), m)
	V = ntuple(x -> zeros(n, n), m)

	for i in eachindex(vcobject)
		copyto!(V[i], vcobject[i].B)
		copyto!(Σ[i], vcobject[i].A)
	end
	return(Σ, V)
end

# for a single Variance Component
# algorithm that will transform z ~ N(0,1)
function SimulateMVN!(z, vc::VarianceComponent)
	#for the ith variance component (VC)
	cholA = vc.CholA # grab (not calculate) the stored Cholesky decomposition of n_traits by n_traits variance component matrix
	cholB = vc.CholB # grab (not calculate) the stored Cholesky decomposition of n_people by n_people variance component matrix

	#Generating MN(0, Sigma)
	# first generate from standard normal
	randn!(z)

	# we want to solve u then v to get the first variance component, v.
	# first matrix vector multiplication using the cholesky decomposed CholA, CholB above
	lmul!(cholB.U, z)
	rmul!(z, cholA.L)

	#second matrix vector multiplication using the he cholesky decomposed CholA, CholB above
	rmul!(z, cholA.U)
	lmul!(cholB.L, z) #multiply on left and save to simulated_trait

	#add the effects of each variance component
	return(z)
end

#this function will call the SimulateMVN! so that I write over z (reuse memory allocation) for potentially many simulations
function SimulateMVN(n_people, n_traits, vc::VarianceComponent)
	#preallocate memory for the returned dataframe simulated_trait once
	z = Matrix{Float64}(undef, n_people, n_traits)

	# calls the function to apply the cholesky decomposition (we have stored in vc object)
	# and transforms/updates z ~ MN(0, vc)
	SimulateMVN!(z, vc::VarianceComponent)
	#returns the allocated and now transformed z
	return(z)
end

##Single Variance Component object
#without computing mean from dataframe and formulas i.e given an evaluated matrix of means
function LMM_trait_simulation(mu, vc::VarianceComponent)
	n_people = size(mu)[1]
	n_traits = size(mu)[2]

	z = SimulateMVN(n_people, n_traits, vc)

	#for each trait add the mean --> MN(mu, Sigma)
	z += mu

	out = DataFrame(z)

	out = names!(out, [Symbol("trait$i") for i in 1:n_traits])

	return out
end

# we note the exclamation is to indicate this function will mutate or override the values that its given
# here z, simulated_trait will get updated but vc will not be touched.
## AGGREGATE MULTIPLE VARIANCE COMPONENTS IN LMM TRAIT OBJECT to creat overall variance
function Aggregate_VarianceComponents!(z, total_variance, vc::Vector{VarianceComponent})
	for i in 1:length(vc)
		SimulateMVN!(z, vc[i])
		#add the effects of each variance component
		total_variance += z
	end

	return total_variance
end

#multiple LMM traits

#without computing mean from dataframe and formulas i.e given an evaluated matrix of means
function LMM_trait_simulation(mu, vc::Vector{VarianceComponent})
	n_people = size(mu)[1]
	n_traits = size(mu)[2]

	#preallocate memory for the returned dataframe simulated_trait
	simulated_trait = zeros(n_people, n_traits)
	z = Matrix{Float64}(undef, n_people, n_traits)

	#using the function above to write over the allocated z and simulated_trait
	Aggregate_VarianceComponents!(z, simulated_trait, vc)
	#for each trait add the mean --> MN(mu, Sigma)
	simulated_trait += mu

	out = DataFrame(simulated_trait)

	out = names!(out, [Symbol("trait$i") for i in 1:n_traits])

	return out
end

max_stars_repo_name: sarah-ji/Tutorials

max_issues_repo_name: sarah-ji/Tutorials

max_forks_repo_name: sarah-ji/Tutorials


text:
module Light

    include("pulse_types.jl")
    include("sech.jl")
    include("OOK.jl")

end


text:
#=
 Copyright 2013-2015 Marco Nehmeier (nehmeier@informatik.uni-wuerzburg.de)
 Copyright 2015-2016 Oliver Heimlich

 Original author: Marco Nehmeier (unit tests in libieeep1788,
                  original license: Apache License 2.0)
 Converted into portable ITL format by Oliver Heimlich with minor corrections.

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation; either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program; if not, see <http://www.gnu.org/licenses/>.
=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
# setrounding(Interval, :correct)
@format full

@testset "minimal_sqr_rev_test" begin
    @test sqr_rev(∅, -∞..∞)[2] == ∅
    @test sqr_rev(Interval(-10.0, -1.0), -∞..∞)[2] == ∅
    @test sqr_rev(Interval(0.0, Inf), -∞..∞)[2] == entireinterval(Float64)
    @test sqr_rev(Interval(0.0, 1.0), -∞..∞)[2] == Interval(-1.0, 1.0)
    @test sqr_rev(Interval(-0.5, 1.0), -∞..∞)[2] == Interval(-1.0, 1.0)
    @test sqr_rev(Interval(-1000.0, 1.0), -∞..∞)[2] == Interval(-1.0, 1.0)
    @test sqr_rev(Interval(0.0, 25.0), -∞..∞)[2] == Interval(-5.0, 5.0)
    @test sqr_rev(Interval(-1.0, 25.0), -∞..∞)[2] == Interval(-5.0, 5.0)
    @test sqr_rev(Interval(0x1.47ae147ae147bp-7, 0x1.47ae147ae147cp-7), -∞..∞)[2] == Interval(-0x1.999999999999bp-4, 0x1.999999999999bp-4)
    @test sqr_rev(Interval(0.0, 0x1.fffffffffffe1p+1), -∞..∞)[2] == Interval(-0x1.ffffffffffff1p+0, 0x1.ffffffffffff1p+0)
end

@testset "minimal_sqr_rev_bin_test" begin
    @test sqr_rev(∅, Interval(-5.0, 1.0))[2] == ∅
    @test sqr_rev(Interval(-10.0, -1.0), Interval(-5.0, 1.0))[2] == ∅
    @test sqr_rev(Interval(0.0, Inf), Interval(-5.0, 1.0))[2] == Interval(-5.0, 1.0)
    @test sqr_rev(Interval(0.0, 1.0), Interval(-0.1, 1.0))[2] == Interval(-0.1, 1.0)
    @test sqr_rev(Interval(-0.5, 1.0), Interval(-0.1, 1.0))[2] == Interval(-0.1, 1.0)
    @test sqr_rev(Interval(-1000.0, 1.0), Interval(-0.1, 1.0))[2] == Interval(-0.1, 1.0)
    @test sqr_rev(Interval(0.0, 25.0), Interval(-4.1, 6.0))[2] == Interval(-4.1, 5.0)
    @test sqr_rev(Interval(-1.0, 25.0), Interval(-4.1, 7.0))[2] == Interval(-4.1, 5.0)
    @test sqr_rev(Interval(1.0, 25.0), Interval(0.0, 7.0))[2] == Interval(1.0, 5.0)
    @test sqr_rev(Interval(0x1.47ae147ae147bp-7, 0x1.47ae147ae147cp-7), Interval(-0.1, Inf))[2] == Interval(-0.1, 0x1.999999999999bp-4)
    @test sqr_rev(Interval(0.0, 0x1.fffffffffffe1p+1), Interval(-0.1, Inf))[2] == Interval(-0.1, 0x1.ffffffffffff1p+0)
end

@testset "minimal_sqr_rev_dec_test" begin
end

@testset "minimal_sqr_rev_dec_bin_test" begin
end

@testset "minimal_abs_rev_test" begin
    @test abs_rev(∅, -∞..∞)[2] == ∅
    @test abs_rev(Interval(-1.1, -0.4), -∞..∞)[2] == ∅
    @test abs_rev(Interval(0.0, Inf), -∞..∞)[2] == entireinterval(Float64)
    @test abs_rev(Interval(1.1, 2.1), -∞..∞)[2] == Interval(-2.1, 2.1)
    @test abs_rev(Interval(-1.1, 2.0), -∞..∞)[2] == Interval(-2.0, 2.0)
    @test abs_rev(Interval(-1.1, 0.0), -∞..∞)[2] == Interval(0.0, 0.0)
    @test abs_rev(Interval(-1.9, 0.2), -∞..∞)[2] == Interval(-0.2, 0.2)
    @test abs_rev(Interval(0.0, 0.2), -∞..∞)[2] == Interval(-0.2, 0.2)
    @test abs_rev(Interval(-1.5, Inf), -∞..∞)[2] == entireinterval(Float64)
end

@testset "minimal_abs_rev_bin_test" begin
    @test abs_rev(∅, Interval(-1.1, 5.0))[2] == ∅
    @test abs_rev(Interval(-1.1, -0.4), Interval(-1.1, 5.0))[2] == ∅
    @test abs_rev(Interval(0.0, Inf), Interval(-1.1, 5.0))[2] == Interval(-1.1, 5.0)
    @test abs_rev(Interval(1.1, 2.1), Interval(-1.0, 5.0))[2] == Interval(1.1, 2.1)
    @test abs_rev(Interval(-1.1, 2.0), Interval(-1.1, 5.0))[2] == Interval(-1.1, 2.0)
    @test abs_rev(Interval(-1.1, 0.0), Interval(-1.1, 5.0))[2] == Interval(0.0, 0.0)
    @test abs_rev(Interval(-1.9, 0.2), Interval(-1.1, 5.0))[2] == Interval(-0.2, 0.2)
end

@testset "minimal_abs_rev_dec_test" begin
end

@testset "minimal_abs_rev_dec_bin_test" begin
end

@testset "minimal_pown_rev_test" begin

end

@testset "minimal_pown_rev_bin_test" begin

end

@testset "minimal_pown_rev_dec_test" begin

end

@testset "minimal_pown_rev_dec_bin_test" begin

end

@testset "minimal_sin_rev_test" begin
    @test sin_rev(∅, -∞..∞)[2] == ∅
    @test sin_rev(Interval(-2.0, -1.1), -∞..∞)[2] == ∅
    @test sin_rev(Interval(1.1, 2.0), -∞..∞)[2] == ∅
    @test sin_rev(Interval(-1.0, 1.0), -∞..∞)[2] == entireinterval(Float64)
    @test sin_rev(Interval(0.0, 0.0), -∞..∞)[2] == entireinterval(Float64)
    @test sin_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), -∞..∞)[2] == entireinterval(Float64)
end

@testset "minimal_sin_rev_bin_test" begin
    @test sin_rev(∅, Interval(-1.2, 12.1))[2] == ∅
    @test sin_rev(Interval(-2.0, -1.1), Interval(-5.0, 5.0))[2] == ∅
    @test sin_rev(Interval(1.1, 2.0), Interval(-5.0, 5.0))[2] == ∅
    @test sin_rev(Interval(-1.0, 1.0), Interval(-1.2, 12.1))[2] == Interval(-1.2, 12.1)
    @test sin_rev(Interval(0.0, 0.0), Interval(-1.0, 1.0))[2] == Interval(0.0, 0.0)
    @test sin_rev(Interval(-0.0, -0.0), Interval(2.0, 2.5))[2] == ∅
    @test sin_rev(Interval(-0.0, -0.0), Interval(3.0, 3.5))[2] == Interval(0x1.921fb54442d18p+1, 0x1.921fb54442d19p+1)
    @test sin_rev(Interval(0x1.fffffffffffffp-1, 0x1p+0), Interval(1.57, 1.58))[2] == Interval(0x1.921fb50442d18p+0, 0x1.921fb58442d1ap+0)
    @test sin_rev(Interval(0.0, 0x1p+0), Interval(-0.1, 1.58))[2] == Interval(0.0, 1.58)
    @test sin_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), Interval(3.14, 3.15))[2] == Interval(0x1.921fb54442d17p+1, 0x1.921fb54442d19p+1)
    @test sin_rev(Interval(-0x1.72cece675d1fdp-52, -0x1.72cece675d1fcp-52), Interval(3.14, 3.15))[2] == Interval(0x1.921fb54442d18p+1, 0x1.921fb54442d1ap+1)
    @test sin_rev(Interval(-0x1.72cece675d1fdp-52, 0x1.1a62633145c07p-53), Interval(3.14, 3.15))[2] == Interval(0x1.921fb54442d17p+1, 0x1.921fb54442d1ap+1)
    @test sin_rev(Interval(0.0, 1.0), Interval(-0.1, 3.15))[2] == Interval(0.0, 0x1.921fb54442d19p+1)
    @test sin_rev(Interval(0.0, 1.0), Interval(-0.1, 3.15))[2] == Interval(-0.0, 0x1.921fb54442d19p+1)
    @test sin_rev(Interval(-0x1.72cece675d1fdp-52, 1.0), Interval(-0.1, 3.15))[2] == Interval(-0x1.72cece675d1fep-52, 0x1.921fb54442d1ap+1)
    @test sin_rev(Interval(-0x1.72cece675d1fdp-52, 1.0), Interval(0.0, 3.15))[2] == Interval(0.0, 0x1.921fb54442d1ap+1)
    @test sin_rev(Interval(0x1.1a62633145c06p-53, 0x1p+0), Interval(3.14, 3.15))[2] == Interval(3.14, 0x1.921fb54442d19p+1)
    @test sin_rev(Interval(-0x1.72cece675d1fdp-52, 0x1p+0), Interval(1.57, 3.15))[2] == Interval(1.57, 0x1.921fb54442d1ap+1)
    @test sin_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), Interval(-Inf, 3.15))[2] == Interval(-Inf, 0x1.921fb54442d19p+1)
    @test sin_rev(Interval(-0x1.72cece675d1fdp-52, -0x1.72cece675d1fcp-52), Interval(3.14, Inf))[2] == Interval(0x1.921fb54442d18p+1, Inf)
end

@testset "minimal_sin_rev_dec_test" begin
end

@testset "minimal_sin_rev_dec_bin_test" begin
end

@testset "minimal_cos_rev_test" begin
    @test cos_rev(∅, -∞..∞)[2] == ∅
    @test cos_rev(Interval(-2.0, -1.1), -∞..∞)[2] == ∅
    @test cos_rev(Interval(1.1, 2.0), -∞..∞)[2] == ∅
    @test cos_rev(Interval(-1.0, 1.0), -∞..∞)[2] == entireinterval(Float64)
    @test cos_rev(Interval(0.0, 0.0), -∞..∞)[2] == entireinterval(Float64)
    @test cos_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), -∞..∞)[2] == entireinterval(Float64)
end

@testset "minimal_cos_rev_bin_test" begin
    @test cos_rev(∅, Interval(-1.2, 12.1))[2] == ∅
    @test cos_rev(Interval(-2.0, -1.1), Interval(-5.0, 5.0))[2] == ∅
    @test cos_rev(Interval(1.1, 2.0), Interval(-5.0, 5.0))[2] == ∅
    @test cos_rev(Interval(-1.0, 1.0), Interval(-1.2, 12.1))[2] == Interval(-1.2, 12.1)
    @test cos_rev(Interval(1.0, 1.0), Interval(-0.1, 0.1))[2] == Interval(0.0, 0.0)
    @test cos_rev(Interval(-1.0, -1.0), Interval(3.14, 3.15))[2] == Interval(0x1.921fb54442d18p+1, 0x1.921fb54442d1ap+1)
    @test cos_rev(Interval(0x1.1a62633145c06p-54, 0x1.1a62633145c07p-54), Interval(1.57, 1.58))[2] == Interval(0x1.921fb54442d17p+0, 0x1.921fb54442d19p+0)
    @test cos_rev(Interval(-0x1.72cece675d1fdp-53, -0x1.72cece675d1fcp-53), Interval(1.57, 1.58))[2] == Interval(0x1.921fb54442d18p+0, 0x1.921fb54442d1ap+0)
    @test cos_rev(Interval(-0x1.72cece675d1fdp-53, 0x1.1a62633145c07p-54), Interval(1.57, 1.58))[2] == Interval(0x1.921fb54442d17p+0, 0x1.921fb54442d1ap+0)
    @test cos_rev(Interval(0x1.1a62633145c06p-54, 1.0), Interval(-2.0, 2.0))[2] == Interval(-0x1.921fb54442d19p+0, 0x1.921fb54442d19p+0)
    @test cos_rev(Interval(0x1.1a62633145c06p-54, 1.0), Interval(0.0, 2.0))[2] == Interval(0.0, 0x1.921fb54442d19p+0)
    @test cos_rev(Interval(-0x1.72cece675d1fdp-53, 1.0), Interval(-0.1, 1.5708))[2] == Interval(-0.1, 0x1.921fb54442d1ap+0)
    @test cos_rev(Interval(-0x1p+0, -0x1.fffffffffffffp-1), Interval(3.14, 3.15))[2] == Interval(0x1.921fb52442d18p+1, 0x1.921fb56442d1ap+1)
    @test cos_rev(Interval(-0x1p+0, -0x1.fffffffffffffp-1), Interval(-3.15, -3.14))[2] == Interval(-0x1.921fb56442d1ap+1, -0x1.921fb52442d18p+1)
    @test cos_rev(Interval(-0x1p+0, -0x1.fffffffffffffp-1), Interval(9.42, 9.45))[2] == Interval(0x1.2d97c7eb321d2p+3, 0x1.2d97c7fb321d3p+3)
    @test cos_rev(Interval(0x1.87996529f9d92p-1, 1.0), Interval(-1.0, 0.1))[2] == Interval(-0x1.6666666666667p-1, 0.1)
    @test cos_rev(Interval(-0x1.aa22657537205p-2, 0x1.14a280fb5068cp-1), Interval(0.0, 2.1))[2] == Interval(0x1.fffffffffffffp-1, 0x1.0000000000001p+1)
    @test cos_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), Interval(-Inf, 1.58))[2] == Interval(-Inf, 0x1.921fb54442d18p+0)
    @test cos_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), Interval(-Inf, 1.5))[2] == Interval(-Inf, -0x1.921fb54442d17p+0)
    @test cos_rev(Interval(-0x1.72cece675d1fdp-52, -0x1.72cece675d1fcp-52), Interval(-1.58, Inf))[2] == Interval(-0x1.921fb54442d1ap+0, Inf)
    @test cos_rev(Interval(-0x1.72cece675d1fdp-52, -0x1.72cece675d1fcp-52), Interval(-1.5, Inf))[2] == Interval(0x1.921fb54442d19p+0, Inf)
end

@testset "minimal_cos_rev_dec_test" begin
end

@testset "minimal_cos_rev_dec_bin_test" begin
end

@testset "minimal_tan_rev_test" begin
    @test tan_rev(∅, -∞..∞)[2] == ∅
    @test tan_rev(Interval(-1.0, 1.0), -∞..∞)[2] == entireinterval(Float64)
    @test tan_rev(Interval(-156.0, -12.0), -∞..∞)[2] == entireinterval(Float64)
    @test tan_rev(Interval(0.0, 0.0), -∞..∞)[2] == entireinterval(Float64)
    @test tan_rev(Interval(0x1.1a62633145c06p-53, 0x1.1a62633145c07p-53), -∞..∞)[2] == entireinterval(Float64)
end

@testset "minimal_tan_rev_bin_test" begin
    @test tan_rev(∅, Interval(-1.5708, 1.5708))[2] == ∅
    @test tan_rev(entireinterval(Float64), Interval(-1.5708, 1.5708))[2] == Interval(-1.5708, 1.5708)
    @test tan_rev(Interval(0.0, 0.0), Interval(-1.5708, 1.5708))[2] == Interval(0.0, 0.0)
    @test tan_rev(Interval(0x1.d02967c31cdb4p+53, 0x1.d02967c31cdb5p+53), Interval(-1.5708, 1.5708))[2] == Interval(-0x1.921fb54442d1bp+0, 0x1.921fb54442d19p+0)
    @test tan_rev(Interval(-0x1.1a62633145c07p-53, -0x1.1a62633145c06p-53), Interval(3.14, 3.15))[2] == Interval(0x1.921fb54442d17p+1, 0x1.921fb54442d19p+1)
    @test tan_rev(Interval(0x1.72cece675d1fcp-52, 0x1.72cece675d1fdp-52), Interval(-3.15, 3.15))[2] == Interval(-0x1.921fb54442d19p+1, 0x1.921fb54442d1ap+1)
    @test tan_rev(Interval(-0x1.d02967c31p+53, 0x1.d02967c31p+53), Interval(-Inf, 1.5707965))[2] == Interval(-Inf, +0x1.921fb82c2bd7fp0)
    @test tan_rev(Interval(-0x1.d02967c31p+53, 0x1.d02967c31p+53), Interval(-1.5707965, Inf))[2] == Interval(-0x1.921fb82c2bd7fp0, Inf)
    @test tan_rev(Interval(-0x1.d02967c31p+53, 0x1.d02967c31p+53), Interval(-1.5707965, 1.5707965))[2] == Interval(-0x1.921fb82c2bd7fp0, +0x1.921fb82c2bd7fp0)
    @test tan_rev(Interval(-0x1.d02967c31cdb5p+53, 0x1.d02967c31cdb5p+53), Interval(-1.5707965, 1.5707965))[2] == Interval(-1.5707965, 1.5707965)
end

@testset "minimal_tan_rev_dec_test" begin
end

@testset "minimal_tan_rev_dec_bin_test" begin
end

@testset "minimal_cosh_rev_test" begin

end

@testset "minimal_cosh_rev_bin_test" begin

end

@testset "minimal_cosh_rev_dec_test" begin

end

@testset "minimal_cosh_rev_dec_bin_test" begin

end
# FactCheck.exitstatus()


text:
@testset "$TEST $G" begin

n = 5
g = CompleteGraph(n, G)
@test count_spanning_trees(g) == 125

g = CompleteGraph(3, G)
@test count_spanning_trees(g) == 3

for n=5:10
    n = 5
    g = PathGraph(n, G)
    @test count_spanning_trees(g) == 1
    add_edge!(g, 1, n)
    @test count_spanning_trees(g) == n
end

g = StarGraph(n, G)
@test count_spanning_trees(g) == 1

g = blockdiag(g, g)
#accept only connected graphs
@test_throws AssertionError count_spanning_trees(g)

# count the empty tree as a spanning tree
g = G()
@test count_spanning_trees(g) == 1
g = G(1)
@test count_spanning_trees(g) == 1

end # testset


text:
#!/usr/bin/env julia

using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, Pose2D, PoseStamped, Vector3
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan

@rosimport gazebo_msgs.msg: ModelState
@rosimport gazebo_msgs.srv: SetModelState, GetModelState, GetWorldProperties

rostypegen()
using geometry_msgs.msg
using std_srvs.srv
using nav_msgs.srv.GetPlan
using gazebo_msgs.msg
using gazebo_msgs.srv

function loop(set_state,get_state)
    loop_rate = Rate(5.0)
    modelName = "robot"

    while !is_shutdown()

        # Get the current position of the Gazebo model
        gs = GetModelStateRequest()
        gs.model_name = modelName
        gs_r = get_state(gs)

        if !gs_r.success
            error(string(" calling /gazebo/get_model_state service: ", gs_r.status_message))
        end

        # Define position to move robot
        vehPose = Pose()
        vehPose.position.x = gs_r.pose.position.x + 0.1

        # Define the robot state
        ms = ModelState()
        ms.model_name = modelName
        ms.pose = vehPose

        # Set the state of the Gazebo model
        ss = SetModelStateRequest()
        ss.model_state = ms
        println("Calling 'gazebo/set_model_state' service...")
        ss_r = set_state(ss)

        if !ss_r.success
            error(string(" calling /gazebo/set_model_state service: ", ss_r.status_message))
        end

        rossleep(loop_rate)
    end
end

function main()
    init_node("rosjl_ex")

    # Set up service to set Gazebo model state
    const set_state = ServiceProxy("/gazebo/set_model_state", SetModelState)
    println("Waiting for 'gazebo/set_model_state' service...")
    wait_for_service("gazebo/set_model_state")

    # Set up service to get Gazebo model state
    const get_state = ServiceProxy("/gazebo/get_model_state",GetModelState)
    println("Waiting for 'gazebo/get_model_state' service...")
    wait_for_service("gazebo/get_model_state")

    loop(set_state,get_state)
end

if ! isinteractive()
    main()
end

#= TODO
# wait for all services to be in some ready state, seehttp://docs.ros.org/electric/api/srs_user_tests/html/prepare__robot__manip__sim_8py_source.html
# start laucnh file in a paused state and have some sin_init file start it

 sim = get_param("/use_sim_time")
 if sim
      loginfo('Waiting until simulated robot is prepared for the task...')

      const get_world = ServiceProxy("/gazebo/get_world_properties",GetWorldProperties)
      println("Waiting for '/gazebo/get_world_properties' service...")
      wait_for_service("/gazebo/get_world_properties")

      wp_r=GetWorldPropertiesRequest()

      wp=get_world(wp_r)

wait_for_message('/sim_robot_init',EmptyMsg)
wait_for_service("/gazebo/pause_physics")

loginfo("Pausing physics")

=#
# world files with obstacles
# LiDAR data processing algorithm
# figure out elipses


text:
struct InvLogitElement{T<:Real} <: Real
    data::T
    logp::T
    nlogitp::T
end
struct InvLogitVec{W,T<:Real} <: AbstractStructVec{W,T}
    data::Vec{W,T}
    logp::Vec{W,T}
    nlogitp::Vec{W,T}
end
const InvLogitValue{T} = Union{InvLogitElement{T},InvLogitVec{<:Any,T}}
@inline invlogitwrap(p::T, logp::T, invlogit::T) where {T<:Number} = InvLogitElement(p, logp, invlogit)
@inline invlogitwrap(p::Vec{W,T}, logp::Vec{W,T}, invlogit::Vec{W,T}) where {W,T<:Number} = InvLogitVec(p, logp, invlogit)
@inline function invlogitwrap(p::AbstractStructVec{W,T}, logp::AbstractStructVec{W,T}, invlogit::AbstractStructVec{W,T}) where {W,T<:Number}
    InvLogitVec(extract_data(p), extract_data(logp), extract_data(invlogit))
end
struct InvLogitStridedPointer{L,T,P<:AbstractStridedPointer{T}} <: AbstractStridedPointer{T}
    ptr::P
    logitptr::Ptr{T}
    @inline InvLogitStridedPointer{L}(ptr::P, logitptr::Ptr{T}) where {L,T,P<:AbstractStridedPointer} = InvLogitStridedPointer{L,T,P}(ptr, logitptr)
end
@inline InvLogitStridedPointer(ptr::P, logitptr::Ptr{T}, ::Val{L}) where {L,T,P<:AbstractStridedPointer} = InvLogitStridedPointer{L,T,P}(ptr, logitptr)
@inline Base.convert(::Type{T1}, a::InvLogitElement{T2}) where {T1,T2} = convert(T1, a.data)
@inline Base.promote_type(::Type{InvLogitElement{T1}}, ::Type{T2}) where {T1,T2} = promote_type(T1, T2)
@inline Base.promote_type(::Type{T2}, ::Type{InvLogitElement{T1}}) where {T1,T2} = promote_type(T1, T2)
@inline Base.promote_type(::Type{InvLogitElement{T1}}, ::Type{InvLogitElement{T2}}) where {T1,T2} = promote_type(T1, T2)
@inline Base.pointer(A::InvLogitStridedPointer) = A.ptr.ptr
@inline VectorizationBase.offset(p::InvLogitStridedPointer, i::Tuple) = offset(p.ptr, i)
@inline function VectorizationBase.vload(ilp::InvLogitStridedPointer{L,T}, i) where {L,T}
    o = offset(ilp, i)
    p = vload(ilp.ptr, o)
    logp = vload(ilp.ptr, o + align(L,T))
    logitp = vload(ilp.logitptr, o)
    invlogitwrap(p, logp, logitp)
end
@inline Base.log(il::InvLogitValue) = il.logp
@inline SLEEFPirates.log1m(il::InvLogitValue) = vadd(il.logp, il.nlogitp)
@inline SLEEFPirates.logit(il::InvLogitValue) = vsub(il.nlogitp)
@inline SLEEFPirates.nlogit(il::InvLogitValue) = il.nlogitp

struct InvLogitArray{S,T,N,X,SN,XN} <: AbstractStrideArray{S,T,N,X,SN,XN,true}
    ptr::Ptr{T}
    logitptr::Ptr{T}
    size::NTuple{SN,UInt32}
    stride::NTuple{XN,UInt32}
end
@generated function InvLogitArray{S}(ptr::Ptr{T}, logitptr::Ptr{T}) where {S,T}
    N, X, L = calc_NXL(S.parameters, T, (S.parameters[1])::Int)
    L = VectorizationBase.align(L, T)
    Expr(:block, Expr(:meta,:inline), :(InvLogitArray{$S,$T,$N,$X,0,0}(ptr, logitptr, tuple(), tuple())))
end
@inline Base.pointer(A::InvLogitArray) = A.ptr
@inline function VectorizationBase.stridedpointer(A::InvLogitArray{S,T,N,X,SN,XN}) where {S,T,N,X,SN,XN}
    InvLogitStridedPointer(stridedpointer(PtrArray(A)), A.logitptr, memory_length_val(A))
end

@inline function constrain(θ::Ptr{Float64}, ::RealScalar{0.0,1.0})
    y = vload(θ)
    p = ninvlogit(y)
    logp = log(p)
    t = Base.FastMath.add_fast(Base.FastMath.mul_fast(2, logp), y)
    t, InvLogitElement(p, logp, y)
end
@inline function constrain_pullback(∇::Ptr{Float64}, θ::Ptr{Float64}, ::RealScalar{0.0,1.0})
    y = vload(θ)
    p = ninvlogit(y)
    logp = log(p)
    t = Base.FastMath.add_fast(Base.FastMath.mul_fast(2, logp), y)
    ile = InvLogitElement(p, logp, y)
    (t, ile), (Zero(), ile)
end
@inline function constrain_double_bounded_array!(c::AbstractStrideMatrix{N,2}, nil::AbstractVector, ::RealArray{S,0.0,1.0,0}) where {N,S}
    t = vzero()
    @avx for i ∈ eachindex(nil)
        y = nil[i]
        p = 1 / (1 + exp(y))
        logp = log(p)
        t += 2*logp + y
        c[i,1] = p
        c[i,2] = logp
    end
    t
end
@inline function stack_pointer_call(::typeof(constrain), sp::StackPointer, θ::Ptr{Float64}, ::RealArray{S,0.0,1.0,0}) where {S}
    sp, e = NoPadFlatPtrViewMulti(sp, S, Val{2}())
    uv = flatvector(NoPadPtrView{S}(θ))
    t = constrain_double_bounded_array!(flatvector(e), uv, RealArray{S,L,U,0}())
    ila = InvLogitArray{S}(pointer(e), pointer(θ))
    sp, (t, ila)
end
@inline function stack_pointer_call(::typeof(constrain_pullback!), sp::StackPointer, ∇::Ptr{Float64}, θ::Ptr{Float64}, ::RealArray{S,0.0,1.0,0}) where {S}
    sp, e = NoPadFlatPtrViewMulti(sp, S, Val{2}())
    uv = flatvector(NoPadPtrView{S}(θ))
    t = constrain_double_bounded_array!(flatvector(e), uv, RealArray{S,L,U,0}())
    ila = InvLogitArray{S}(pointer(e), pointer(θ))
    sp, ((t, ila), (NoPadPtrView{S}(∇), ila))
end


text:
@testset "I473" begin
    s = """
    r = [3,2,1]
    r[r]
    """
    msgs = lintstr(s)
    @test messageset(msgs) == Set([:I473])
    @test occursin("value at position #1 is the referenced r", msgs[1].message)
    @test occursin("OK if it represents permutations", msgs[1].message)
end

@testset "E434" begin
    msgs = lintstr("""
    r = [3,2,1]
    r[1;r]
    """)
    @test messageset(msgs) == Set([:E434])
    @test occursin("value at position #2 is the referenced r. Possible typo?", msgs[1].message)
end


text:
function AxisIndices.axes_type(::Type{<:NamedDimsArray{L,T,N,A}}) where {L,T,N,A}
    return AxisIndices.axes_type(A)
end

"""
    NIArray()
"""
const NIArray{names,T,N,P<:AbstractArray{T,N},Ax} = NamedDimsArray{names,T,N,AxisIndicesArray{T,N,P,Ax}}
NIArray(x::AbstractArray; kwaxes...) = NIArray(x, kwaxes.data)
function NIArray(x::AbstractArray, namedaxes::NamedTuple{L}) where {L}
    return NamedDimsArray{L}(AxisIndicesArray(x, values(namedaxes)))
end

"""
    IMArray()
"""
const IMArray{T,N,P<:AbstractArray{T,N},Ax,M} = AxisIndicesArray{T,N,ImageMeta{T,N,P,M},Ax}
IMArray(x::AbstractArray, args...; metadata=Metadata()) = IMArray(x, Tuple(args), metadata)
IMArray(x::AbstractArray, axs::Tuple, metadata) = AxisIndicesArray(ImageMeta(x, metadata), axs)

"""
    MetaArray

"""
const MetaArray{T,N,P<:AbstractArray{T,N},M<:AbstractMetadata} = ImageMeta{T,N,P,M}
MetaArray(x::AbstractArray, metadata=Metadata()) = ImageMeta(x, metadata)

"""
    NIMArray(x, args...; metadata, kwargs...)

Returns an array with named dimensions, an `AbstractAxis` for each indice, and metadata)
"""
const NIMArray{names,T,N,P<:AbstractArray{T,N},Ax,M} = NamedDimsArray{names,T,N,AxisIndicesArray{T,N,ImageMeta{T,N,P,M},Ax}}

function NIMArray(x::AbstractArray; metadata=Metadata(), kwaxes...)
    return NIMArray(x, kwaxes.data, metadata)
end

function NIMArray(x::AbstractArray, namedaxes::NamedTuple{L}, metadata) where {L}
    return NamedDimsArray{L}(IMArray(x, values(namedaxes), metadata))
end

function NIMArray(x::AbstractArray{T,N}, args...; metadata::AbstractDict{Symbol,Any}=Metadata(), kwargs...) where {T,N}
    if isempty(args)
        if isempty(kwargs)
            return _NIMArray(x, metadata)
        else
            return _NIMArray(x, values(kwargs), metadata)
        end
    elseif isempty(kwargs)
        return _NIMArray(x, args, metadata)
    else
        error("Indices can only be specified by keywords or additional arguments after the parent array, not both.")
    end
end

function _NIMArray(x::AbstractArray{T,N}, namedindices::NamedTuple{L}, metadata) where {T,N,L}
    return NamedDimsArray{L}(AxisIndicesArray(ImageMeta(x, metadata), values(namedindices)))
end

function _NIMArray(x::AbstractArray{T,N}, indices::Tuple, metadata) where {T,N,L}
    return NamedDimsArray{ntuple(_ -> :_, Val(N))}(AxisIndicesArray(ImageMeta(x, metadata), indices))
end

function _NIMArray(x::AbstractArray{T,N}, metadata) where {T,N,L}
    return NamedDimsArray{ntuple(_ -> :_, Val(N))}(AxisIndicesArray(ImageMeta(x, metadata)))
end

function _NIMArray(x::AbstractArray{T,N}, axs::Tuple{Vararg{Symbol}}, metadata) where {T,N}
    return NamedDimsArray{axs}(AxisIndicesArray(ImageMeta(x, metadata)))
end

const NamedMappedArray{L,T,N,A,F} = Union{<:ReadonlyMappedArray{T,N,NamedDimsArray{L,T,N,A},F},<:MappedArray{T,N,NamedDimsArray{L,T,N,A},F}}

NamedDims.dimnames(::Type{<:NamedMappedArray{L,T,N,A,F}}) where {L,T,N,A,F} = L
function AxisIndices.axes_type(::Type{<:NamedMappedArray{L,T,N,A,F}}) where {L,T,N,A,F}
    return AxisIndices.axes_type(A)
end

@inline function Base.axes(nma::NamedMappedArray{L,T,N,A,F}, i::Symbol) where {L,T,N,A,F}
    return axes(parent(nma), dim(L, i))
end

# TODO move to DimNames
@inline named_axes(x::X) where {X} = NamedTuple{dimnames(X)}(axes(x))


text:
import Base: +, -, isequal, isless, print, show
import Base.Dates: AbstractTime, Second, value

# Note: The IANA time zone database rounds offset precision to the nearest second
# See "America/New_York" notes in tzdata file "northamerica" for an example.
"""
    UTCOffset

A `UTCOffset` is an amount of time subtracted from or added to UTC to get the current
local time – whether it's standard time or daylight saving time.
"""
immutable UTCOffset <: AbstractTime
    std::Second  # Standard time offset from UTC in seconds
    dst::Second  # Daylight saving time offset in seconds

    function UTCOffset(std_offset::Second, dst_offset::Second=Second(0))
        new(std_offset, dst_offset)
    end
end

function UTCOffset(std_offset::Integer, dst_offset::Integer=0)
    UTCOffset(Second(std_offset), Second(dst_offset))
end

value(offset::UTCOffset) = value(offset.std + offset.dst)

(+)(dt::DateTime, offset::UTCOffset) = dt + (offset.std + offset.dst)
(-)(dt::DateTime, offset::UTCOffset) = dt - (offset.std + offset.dst)

# Determines if the given `UTCOffset` is an offset for daylight saving time.
isdst(offset::UTCOffset) = offset.dst != Second(0)

# Two `UTCOffset`s can be considered equal if the total offset is the same and they are
# both either offsets for standard time or daylight saving time.
function isequal(x::UTCOffset, y::UTCOffset)
    x == y || value(x) == value(y) && isdst(x) == isdst(y)
end
isless(x::UTCOffset, y::UTCOffset) = isless(value(x), value(y))

function offset_string(seconds::Second, iso8601::Bool=false)
    v = value(seconds)
    h, v = divrem(v, 3600)
    m, s  = divrem(abs(v), 60)

    if !iso8601 && m == 0 && s == 0
        return @sprintf("%+02d", h)
    elseif s == 0
        return @sprintf("%+03d:%02d", h, m)
    else
        return @sprintf("%+03d:%02d:%02d", h, m, s)  # Not in the ISO 8601 standard
    end
end
function offset_string(offset::UTCOffset, iso8601::Bool=false)
    offset_string(offset.std + offset.dst, iso8601)
end

print(io::IO, o::UTCOffset) = print(io, offset_string(o, true))
function show(io::IO, o::UTCOffset)
    # Show DST as a separate offset since we want to distinguish between normal hourly
    # daylight saving time offsets and exotic DST offsets (e.g. midsummer time).
    print(io, "UTC", offset_string(o.std), "/", offset_string(o.dst))
end


text:
import Pkg
Pkg.add("Documenter")
using Documenter
using XLStats
push!(LOAD_PATH,"../src/")
makedocs(
    modules=[XLStats],
    sitename="XLStats.jl",
    pages = [
        "Home" => "index.md",
    ]
)
deploydocs(
    repo = "github.com/m3g/XLStats.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#" ],
)


text:
const IntChar = UInt16
const IntSmallString = UInt256

# XXX: remove smallLength as a global constant, bad for type inference
const smallLength = 16
const SmallStringStorage = SVector{smallLength,IntChar}
const MSmallStringStorage = MVector{smallLength,IntChar}

# Similar types are implemented in various packages:
# https://github.com/JuliaString/ShortStrings.jl
# https://github.com/JuliaComputing/FixedSizeStrings.jl
# https://gist.github.com/SimonDanisch/02e74622e0577f199c1b1a8a65390c24#file-fixedstring-jl
# https://github.com/JuliaStrings/StringViews.jl
# https://discourse.julialang.org/t/way-to-make-sharedarray-over-fixed-length-strings/7082
# https://github.com/djsegal/FixedLengthStrings.jl
# TODO: make this more generic by parametrizing over the length and Char size. Also, store the length of the string.
struct SmallString
  data::SmallStringStorage

  SmallString(sv::SmallStringStorage) = new(sv)

  function SmallString()
    store = SmallStringStorage(ntuple(_ -> IntChar(0), Val(smallLength)))
    return new(store)
  end
end

const Tag = SmallString

data(ss::SmallString) = ss.data

Base.eltype(ss::SmallString) = eltype(data(ss))

function SmallString(str)
  length(str) > smallLength &&
    error("String is too long for SmallString. Maximum length is $smallLength.")
  mstore = MSmallStringStorage(ntuple(_ -> IntChar(0), Val(smallLength)))
  for (n, c) in enumerate(str)
    mstore[n] = IntChar(c)
  end
  return SmallString(SmallStringStorage(mstore))
end

SmallString(s::SmallString) = SmallString(data(s))

Base.getindex(s::SmallString, n::Int) = getindex(s.data, n)

function Base.setindex(s::SmallString, val, n::Int)
  return SmallString(StaticArrays.setindex(s.data, val, n))
end

# TODO: rename to `isempty`
isnull(s::SmallString) = @inbounds s[1] == IntChar(0)

function Base.vcat(s1::SmallString, s2::SmallString)
  v = MSmallStringStorage(ntuple(_ -> IntChar(0), Val(smallLength)))
  n = 1
  while n <= smallLength && s1[n] != IntChar(0)
    v[n] = s1[n]
    n += 1
  end
  N1 = n - 1
  n2 = 1
  while n2 <= smallLength && s2[n2] != IntChar(0)
    v[n] = s2[n2]
    n += 1
    n2 += 1
  end
  return SmallString(SmallStringStorage(v))
end

function SmallString(i::IntSmallString)
  mut_is = MVector{1,IntSmallString}(ntoh(i))
  p = convert(Ptr{SmallStringStorage}, pointer_from_objref(mut_is))
  return SmallString(unsafe_load(p))
end

function cast_to_uint(store)
  mut_store = MSmallStringStorage(store)
  storage_begin = convert(Ptr{IntSmallString}, pointer_from_objref(mut_store))
  return ntoh(unsafe_load(storage_begin))
end

function IntSmallString(s::SmallString)
  return cast_to_uint(s.data)
end

function isint(s::SmallString)::Bool
  ndigits = 1
  while ndigits <= smallLength && s[ndigits] != IntChar(0)
    cur_char = Char(s[ndigits])
    !isdigit(cur_char) && return false
    ndigits += 1
  end
  return true
end

# Here we use the StaticArrays comparison
Base.:(==)(s1::SmallString, s2::SmallString) = (s1.data == s2.data)
Base.isless(s1::SmallString, s2::SmallString) = isless(s1.data, s2.data)

########################################################
# Here are alternative SmallString comparison implementations
#

#Base.isless(a::SmallString,b::SmallString) = cast_to_uint(a) < cast_to_uint(b)
#Base.:(==)(a::SmallString,b::SmallString) = cast_to_uint(a) == cast_to_uint(b)

# Here we use the c-function memcmp (used in Julia string comparison):
#function Base.cmp(a::SmallString, b::SmallString)
#  return ccall(:memcmp, Int32, (Ptr{IntChar}, Ptr{IntChar}, IntChar),
#               a.data, b.data, IntChar(8))
#end
#Base.isless(a::SmallString,b::SmallString) = cmp(a, b) < 0
#Base.:(==)(a::SmallString,b::SmallString) = cmp(a, b) == 0

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

function Base.String(s::SmallString)
  n = 1
  while n <= smallLength && s[n] != IntChar(0)
    n += 1
  end
  len = n - 1
  return String(Char.(s.data[1:len]))
end

function Base.show(io::IO, s::SmallString)
  n = 1
  while n <= smallLength && s[n] != IntChar(0)
    print(io, Char(s[n]))
    n += 1
  end
end

function readcpp(io::IO, ::Type{SmallString}; kwargs...)
  format = get(kwargs, :format, "v3")
  s = SmallString()
  if format == "v3"
    for n in 1:7
      c = read(io, Char)
      s = setindex(s, c, n)
    end
  else
    throw(ArgumentError("read SmallString: format=$format not supported"))
  end
  return s
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, "../src/"))

import Documenter, BibTeX

Documenter.makedocs(
    sitename = "BibTeX.jl",
)

Documenter.deploydocs(
    repo = "github.com/JuliaTeX/BibTeX.jl.git",
    target = "build",
    deps = nothing,
    make = nothing
)


text:
#
#  1. TIMESTEP
#

"""
	gettime(ts::FixedTimestep)

Return the time (year) represented by Timestep `ts`
"""
function gettime(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	return FIRST + (ts.t - 1) * STEP
end

"""
	gettime(ts::VariableTimestep)

Return the time (year) represented by Timestep `ts`
"""
function gettime(ts::VariableTimestep)
	return ts.current
end

"""
	is_time(ts::AbstractTimestep, t::Int)

Return true or false, true if the current time (year) for `ts` is `t`
"""
function is_time(ts::AbstractTimestep, t::Int)
	return gettime(ts) == t
end

"""
	is_first(ts::AbstractTimestep)

Return true or false, true if `ts` is the first timestep to be run.
"""
function is_first(ts::AbstractTimestep)
	return ts.t == 1
end

"""
	is_timestep(ts::AbstractTimestep, t::Int)

Return true or false, true if `ts` timestep is step `t`.
"""
function is_timestep(ts::AbstractTimestep, t::Int)
	return ts.t == t
end

"""
	is_last(ts::FixedTimestep)

Return true or false, true if `ts` is the last timestep to be run.
"""
function is_last(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	return gettime(ts) == LAST
end

"""
	is_last(ts::VariableTimestep)

Return true or false, true if `ts` is the last timestep to be run.  Note that you may
run `next_timestep` on `ts`, as ths final timestep has not been run through yet.
"""
function is_last(ts::VariableTimestep{TIMES}) where {TIMES}
	return gettime(ts) == TIMES[end]
end

function finished(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	return gettime(ts) > LAST
end

function finished(ts::VariableTimestep{TIMES}) where {TIMES}
	return gettime(ts) > TIMES[end]
end

function next_timestep(ts::FixedTimestep{FIRST, STEP, LAST}) where {FIRST, STEP, LAST}
	if finished(ts)
			error("Cannot get next timestep, this is last timestep.")
	end
	return FixedTimestep{FIRST, STEP, LAST}(ts.t + 1)
end

function next_timestep(ts::VariableTimestep{TIMES}) where {TIMES}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	end
	return VariableTimestep{TIMES}(ts.t + 1)
end

function Base.:-(ts::FixedTimestep{FIRST, STEP, LAST}, val::Int) where {FIRST, STEP, LAST}
	if is_first(ts)
		error("Cannot get previous timestep, this is first timestep.")
	elseif ts.t - val <= 0
		error("Cannot get requested timestep, precedes first timestep.")
	end
	return FixedTimestep{FIRST, STEP, LAST}(ts.t - val)
end

function Base.:-(ts::VariableTimestep{TIMES}, val::Int) where {TIMES}
	if is_first(ts)
		error("Cannot get previous timestep, this is first timestep.")
	elseif ts.t - val <= 0
		error("Cannot get requested timestep, precedes first timestep.")
	end
	return VariableTimestep{TIMES}(ts.t - val)
end

function Base.:+(ts::FixedTimestep{FIRST, STEP, LAST}, val::Int) where {FIRST, STEP, LAST}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	elseif gettime(ts) + val > LAST + 1
		error("Cannot get requested timestep, exceeds last timestep.")
	end
	new_ts = FixedTimestep{FIRST, STEP, LAST}(ts.t + val)

end

function Base.:+(ts::VariableTimestep{TIMES}, val::Int) where {TIMES}
	if finished(ts)
		error("Cannot get next timestep, this is last timestep.")
	elseif gettime(ts) + val > TIMES[end] + 1
		error("Cannot get requested timestep, exceeds last timestep.")
	end
	new_ts = VariableTimestep{TIMES}(ts.t + val)
end

#
#  2. CLOCK
#

function timestep(c::Clock)
	return c.ts
end

function time_index(c::Clock)
	return c.ts.t
end

"""
	gettime(c::Clock)

Return the current time of the timestep held by the `c` clock.
"""
function gettime(c::Clock)
	return gettime(c.ts)
end

function advance(c::Clock)
	c.ts = next_timestep(c.ts)
	nothing
end

function finished(c::Clock)
	return finished(c.ts)
end

#
# 3.  TimestepVector and TimestepMatrix
#

#
# 3a.  General
#

# Get a timestep array of type T with N dimensions. Time labels will match those from the time dimension in md
function get_timestep_array(md::ModelDef, T, N, value)

	if isuniform(md)
        first, stepsize = first_and_step(md)
        return TimestepArray{FixedTimestep{first, stepsize}, T, N}(value)
    else
        TIMES = (time_labels(md)...,)
        return TimestepArray{VariableTimestep{TIMES}, T, N}(value)
    end
end

const AnyIndex = Union{Int, Vector{Int}, Tuple, Colon, OrdinalRange}

# TBD: can it be reduced to this?
# const AnyIndex = Union{Int, AbstractRange}

#
# 3b. TimestepVector
#

function Base.getindex(v::TimestepVector{FixedTimestep{FIRST, STEP}, T}, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST}
	data = v.data[ts.t]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

function Base.getindex(v::TimestepVector{VariableTimestep{TIMES}, T}, ts::VariableTimestep{TIMES}) where {T, TIMES}
	data = v.data[ts.t]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

function Base.getindex(v::TimestepVector{FixedTimestep{D_FIRST, STEP}, T}, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	data = v.data[t]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

function Base.getindex(v::TimestepVector{VariableTimestep{D_FIRST}, T}, ts::VariableTimestep{T_FIRST}) where {T, D_FIRST, T_FIRST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	data = v.data[t]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

# int indexing version supports old-style components and internal functions, not
# part of the public API

function Base.getindex(v::TimestepVector{FixedTimestep{FIRST, STEP}, T}, i::AnyIndex) where {T, FIRST, STEP}
	return v.data[i]
end

function Base.getindex(v::TimestepVector{VariableTimestep{TIMES}, T}, i::AnyIndex) where {T, TIMES}
	return v.data[i]
end

function Base.setindex!(v::TimestepVector{FixedTimestep{FIRST, STEP}, T}, val, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, FIRST, STEP, LAST}
	setindex!(v.data, val, ts.t)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{TIMES}, T}, val, ts::VariableTimestep{TIMES}) where {T, TIMES}
	setindex!(v.data, val, ts.t)
end

function Base.setindex!(v::TimestepVector{FixedTimestep{D_FIRST, STEP}, T}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	setindex!(v.data, val, t)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{D_FIRST}, T}, val, ts::VariableTimestep{T_FIRST}) where {T, D_FIRST, T_FIRST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	setindex!(v.data, val, t)
end

# int indexing version supports old-style components and internal functions, not
# part of the public API

function Base.setindex!(v::TimestepVector{FixedTimestep{Start, STEP}, T}, val, i::AnyIndex) where {T, Start, STEP}
	setindex!(v.data, val, i)
end

function Base.setindex!(v::TimestepVector{VariableTimestep{TIMES}, T}, val, i::AnyIndex) where {T, TIMES}
	setindex!(v.data, val, i)
end

function Base.length(v::TimestepVector)
	return length(v.data)
end

Base.lastindex(v::TimestepVector) = length(v)

#
# 3c. TimestepMatrix
#

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, ts::FixedTimestep{FIRST, STEP, LAST}, i::AnyIndex) where {T, FIRST, STEP, LAST}
	data = mat.data[ts.t, i]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, ts::VariableTimestep{TIMES}, i::AnyIndex) where {T, TIMES}
	data = mat.data[ts.t, i]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

function Base.getindex(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T}, ts::FixedTimestep{T_FIRST, STEP, LAST}, i::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	data = mat.data[t, i]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{D_FIRST}, T}, ts::VariableTimestep{T_FIRST}, i::AnyIndex) where {T, D_FIRST, T_FIRST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	data = mat.data[t, i]
	if data === missing
		error("Cannot get index; data is missing. You may have tried to access a value that has not yet been computed.")
	else
		return data
	end
end

# int indexing version supports old-style components and internal functions, not
# part of the public API

function Base.getindex(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, idx1::AnyIndex, idx2::AnyIndex) where {T, FIRST, STEP}
	return mat.data[idx1, idx2]
end

function Base.getindex(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, idx1::AnyIndex, idx2::AnyIndex) where {T, TIMES}
	return mat.data[idx1, idx2]
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, val, ts::FixedTimestep{FIRST, STEP, LAST}, idx::AnyIndex) where {T, FIRST, STEP, LAST}
	setindex!(mat.data, val, ts.t, idx)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, val, ts::VariableTimestep{TIMES}, idx::AnyIndex) where {T, TIMES}
	setindex!(mat.data, val, ts.t, idx)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{D_FIRST, STEP}, T}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (T_FIRST - D_FIRST) / STEP)
	setindex!(mat.data, val, t, idx)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{D_FIRST}, T}, val, ts::VariableTimestep{T_FIRST}, idx::AnyIndex) where {T, D_FIRST, T_FIRST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	setindex!(mat.data, val, t, idx)
end

# int indexing version supports old-style components and internal functions, not
# part of the public API

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, val, idx1::Int, idx2::Int) where {T, FIRST, STEP}
	setindex!(mat.data, val, idx1, idx2)
end

function Base.setindex!(mat::TimestepMatrix{FixedTimestep{FIRST, STEP}, T}, val, idx1::AnyIndex, idx2::AnyIndex) where {T, FIRST, STEP}
	mat.data[idx1,idx2] .= val
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, val, idx1::Int, idx2::Int) where {T, TIMES}
	setindex!(mat.data, val, idx1, idx2)
end

function Base.setindex!(mat::TimestepMatrix{VariableTimestep{TIMES}, T}, val, idx1::AnyIndex, idx2::AnyIndex) where {T, TIMES}
	mat.data[idx1,idx2] .= val
end

#
# 4. TimestepArray methods
#

Base.fill!(obj::TimestepArray, value) = fill!(obj.data, value)

Base.size(obj::TimestepArray) = size(obj.data)

Base.size(obj::TimestepArray, i::Int) = size(obj.data, i)

Base.ndims(obj::TimestepArray{T_ts, T, N}) where {T_ts,T, N} = N

Base.eltype(obj::TimestepArray{T_ts, T, N}) where {T_ts,T, N} = T

first_period(obj::TimestepArray{FixedTimestep{FIRST,STEP}, T, N}) where {FIRST, STEP, T, N} = FIRST
first_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N}) where {TIMES, T, N} = TIMES[1]

last_period(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}) where {FIRST, STEP,T, N} = (FIRST + (size(obj, 1) - 1) * STEP)
last_period(obj::TimestepArray{VariableTimestep{TIMES}, T, N}) where {TIMES,T, N} = TIMES[end]

time_labels(obj::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}) where {FIRST, STEP, T, N} = collect(FIRST:STEP:(FIRST + (size(obj, 1) - 1) * STEP))
time_labels(obj::TimestepArray{VariableTimestep{TIMES}, T, N}) where {TIMES, T, N} = collect(TIMES)

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, ts::FixedTimestep{FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, FIRST, STEP, LAST}
	return arr.data[ts.t, idxs...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, ts::VariableTimestep{TIMES}, idxs::AnyIndex...) where {T, N, TIMES}
	return arr.data[ts.t, idxs...]
end

function Base.getindex(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, ts::FixedTimestep{T_FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
	t = Int(ts.t + (FIRST - TIMES[1]) / STEP)
	return arr.data[t, idxs...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, ts::VariableTimestep{T_FIRST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	return arr.data[t, idxs...]
end

# int indexing version supports old-style components and internal functions, not
# part of the public API; first index is Int or Range, rather than a Timestep

function Base.getindex(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, FIRST, STEP}
	return arr.data[idx1, idx2, idxs...]
end

function Base.getindex(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, TIMES}
	return arr.data[idx1, idx2, idxs...]
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, val, ts::FixedTimestep{FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, FIRST, STEP, LAST}
	setindex!(arr.data, val, ts.t, idxs...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, val, ts::VariableTimestep{TIMES}, idxs::AnyIndex...) where {T, N, TIMES}
	setindex!(arr.data, val, ts.t, idxs...)
end

function Base.setindex!(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, val, ts::FixedTimestep{T_FIRST, STEP, LAST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	setindex!(arr.data, val, t, idxs...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, val, ts::VariableTimestep{T_FIRST}, idxs::AnyIndex...) where {T, N, D_FIRST, T_FIRST}
	t = ts.t + findfirst(isequal(T_FIRST[1]), D_FIRST) - 1
	setindex!(arr.data, val, t, idxs...)
end

# int indexing version supports old-style components and internal functions, not
# part of the public API; first index is Int or Range, rather than a Timestep

function Base.setindex!(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, val, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, FIRST, STEP}
	setindex!(arr.data, val, idx1, idx2, idxs...)
end

function Base.setindex!(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, val, idx1::AnyIndex, idx2::AnyIndex, idxs::AnyIndex...) where {T, N, TIMES}
	setindex!(arr.data, val, idx1, idx2, idxs...)
end

"""
	hasvalue(arr::TimestepArray, ts::FixedTimestep)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts`.
"""
function hasvalue(arr::TimestepArray{FixedTimestep{FIRST, STEP}, T, N}, ts::FixedTimestep{FIRST, STEP, LAST}) where {T, N, FIRST, STEP, LAST}
	return 1 <= ts.t <= size(arr, 1)
end

"""
	hasvalue(arr::TimestepArray, ts::VariableTimestep)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts`.
"""
function hasvalue(arr::TimestepArray{VariableTimestep{TIMES}, T, N}, ts::VariableTimestep{TIMES}) where {T, N, TIMES}
	return 1 <= ts.t <= size(arr, 1)
end

function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N}, ts::FixedTimestep{T_FIRST, STEP, LAST}) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
	return D_FIRST <= gettime(ts) <= last_period(arr)
end

function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N}, ts::VariableTimestep{T_FIRST}) where {T, N, T_FIRST, D_FIRST}
	return D_FIRST[1] <= gettime(ts) <= last_period(arr)
end

"""
	hasvalue(arr::TimestepArray, ts::FixedTimestep, idxs::Int...)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts` within
indices `idxs`. Used when Array and Timestep have different FIRST, validating all dimensions.
"""
function hasvalue(arr::TimestepArray{FixedTimestep{D_FIRST, STEP}, T, N},
	ts::FixedTimestep{T_FIRST, STEP, LAST},
	idxs::Int...) where {T, N, D_FIRST, T_FIRST, STEP, LAST}
	return D_FIRST <= gettime(ts) <= last_period(arr) && all([1 <= idx <= size(arr, i) for (i, idx) in enumerate(idxs)])
end

"""
	hasvalue(arr::TimestepArray, ts::VariableTimestep, idxs::Int...)

Return `true` or `false`, `true` if the TimestepArray `arr` contains the Timestep `ts` within
indices `idxs`. Used when Array and Timestep different TIMES, validating all dimensions.
"""
function hasvalue(arr::TimestepArray{VariableTimestep{D_FIRST}, T, N},
	ts::VariableTimestep{T_FIRST},
	idxs::Int...) where {T, N, D_FIRST, T_FIRST}

	return D_FIRST[1] <= gettime(ts) <= last_period(arr) && all([1 <= idx <= size(arr, i) for (i, idx) in enumerate(idxs)])
end


text:
module FiniteFields

##############################################################################
##
## TODO
##
##############################################################################

# - Should remove module and rely on package-level module as per convention
# - Methods e.g., promote, convert, macro +, -, *, /

##############################################################################
##
## Exports
##
##############################################################################

export isIrreducible,
       FiniteField,
       FiniteFieldElement,
       inverse

##############################################################################
##
## Implementation
##
##############################################################################

using Crypto
using Crypto.Integers
using Crypto.Polynomials

import Base.show,
       Base.length,
       Base.convert,
       Base.divrem,
       Base.promote_rule,
       Base.promote,
       Base.abs,
       Base.div,
       Base.rem,
       Base.zero,
       Base.zeros

# TODO: Should make another type, not Number
type FiniteField{P, M}
  polynomial::Polynomial{IntegerMod{P}}
  polynomialModulus::Polynomial{IntegerMod{P}}

  function FiniteField()
    polyMod::Polynomial{IntegerMod{P}} = generateIrreducibleModularPolynomial(P, M)
    new(Polynomial{IntegerMod{P}}(), polyMod)
  end
  function FiniteField(poly::Number)
    polyMod::Polynomial{IntegerMod{P}} = generateIrreducibleModularPolynomial(P, M)
    new(Polynomial{IntegerMod{P}}([poly]), polyMod)
  end
  function FiniteField(poly::Array)
    polyMod::Polynomial{IntegerMod{P}} = generateIrreducibleModularPolynomial(P, M)
    new(Polynomial{IntegerMod{P}}([poly]), polyMod)
  end
  function FiniteField(poly::Polynomial)
    polyMod::Polynomial{IntegerMod{P}} = generateIrreducibleModularPolynomial(P, M)
    new(poly, polyMod)
  end
end

function FiniteField(p::Integer, m::Integer)
  return FiniteField{p, m}
end

# convert{P, M}(::Type{FiniteField{P, M}}, n::Number) = FiniteField{P, M}(n)
# promote_rule{P, M}(::Type{FiniteField{P, M}}, ::Type{Number} ) = FiniteField{P, M}

### Functions ###

function isIrreducible{P}(poly::Polynomial{IntegerMod{P}})
  Zp = IntegerMod{P}
  x = Polynomial([Zp(0), Zp(1)])
  powerTerm = x
  isUnit(pol::Polynomial) = degree(pol) == 0

  for _ in 1:int(degree(poly) / 2)
    powerTerm = powerTerm ^ P
    powerTerm = rem(powerTerm, poly)
    gcdModP = gcd(poly, powerTerm - x)
    if !isUnit(gcdModP)
      return false
    end
  end

  return true
end

function generateIrreducibleModularPolynomial(p::Integer, m::Integer)
  modulus = p
  degree = m

  Zp = IntegerMod{modulus}
  while true
    firstRand::Array{Zp} = [Zp(rand(0:modulus - 1)) for i in 1:degree]
    append!(firstRand, [Zp(1)])

    p = Polynomial{Zp}(firstRand)
    if degree <= 1 || isIrreducible(p)
      return p
    end
  end
end

### Operators on FiniteField and FiniteFieldElement ###

function show{P, M}(io::IO, f::FiniteField{P, M})
  print(io, string(f.polynomial, ' ', '\u2208', ' ', "F_{", P, "^", M, "}"))
end

function +{P, M}(a::FiniteField{P, M}, b::FiniteField{P, M})
  return FiniteField{P, M}(a.polynomial + b.polynomial)
end

function +{P, M}(a::FiniteField{P, M}, b::Number)
  return FiniteField{P, M}(a.polynomial + Polynomial{IntegerMod{P}}(b))
end

function +{P, M}(a::Number, b::FiniteField{P, M})
  return FiniteField{P, M}(Polynomial{IntegerMod{P}(a)} + b.polynomial)
end

function -{P, M}(a::FiniteField{P, M}, b::FiniteField{P, M})
  return FiniteField{P, M}(a.polynomial - b.polynomial)
end

function -{P, M}(a::FiniteField{P, M}, b::Number)
  return FiniteField{P, M}(a.polynomial - Polynomial{IntegerMod{P}}(b))
end

function -{P, M}(a::Number, b::FiniteField{P, M})
  return FiniteField{P, M}(Polynomial{IntegerMod{P}(a)} - b.polynomial)
end

function *{P, M}(a::FiniteField{P, M}, b::FiniteField{P, M})
  return FiniteField{P, M}(a.polynomial * b.polynomial)
end

function *{P, M}(a::FiniteField{P, M}, b::Number)
  return FiniteField{P, M}(a.polynomial * Polynomial{IntegerMod{P}}(b))
end

function *{P, M}(a::Number, b::FiniteField{P, M})
  return FiniteField{P, M}(Polynomial{IntegerMod{P}(a)} * b.polynomial)
end

function =={P, M}(a::FiniteField{P, M}, b::FiniteField{P, M})
  return FiniteField{P, M}(a.polynomial == b.polynomial)
end

function ^{P, M}(a::FiniteField{P, M}, b::Integer)
  return FiniteField{P, M}(a.polynomial ^ b)
end

function -{P, M}(a::FiniteField{P, M})
  return FiniteField{P, M}(-a.polynomial)
end

function divrem{P, M}(a::FiniteField{P, M}, b::FiniteField{P, M})
  q, r = divrem(a.polynomial, b.polynomial)
  return (FiniteField{P, M}(q), FiniteField{P, M}(r))
end

function inverse{P, M}(a::FiniteField{P, M})
  if a.polynomial == Polynomial{IntegerMod{P}}()
    return Inf # TODO: Should this be an error?
  end

  x, y, d = extendedEuclideanAlgorithm(a.polynomial, a.polynomialModulus)
  # ERROR: inverse of d.coefficients?
  return FiniteField{P, M}(x) * FiniteField{P, M}(inverse(d.coefficients[1]))
end

end # module FiniteFields

# type FiniteFieldElement{P, M} <: Number
#   element::Polynomial{IntegerMod{P}}
#   f::FiniteField{P,M}
#   function FiniteFieldElement(e::Polynomial{IntegerMod{P}}, f::FiniteField{P,M})
#     mod = rem(e,f.polynomialModulus)
#     new(mod,f)
#   end
# end

# function FiniteFieldElement{P,M}(e::FiniteFieldElement{P,M}, f::FiniteField{P,M})
#   FiniteFieldElement{P,M}(e.element,f)
# end

# function FiniteFieldElement{P,M}(i::IntegerMod{P}, f::FiniteField{P,M})
#   FiniteFieldElement{P,M}(Polynomial(i),f)
# end

# function FiniteFieldElement{P,M}(i::Int64, f::FiniteField{P,M})
#   FiniteFieldElement{P,M}(Polynomial(makeModular(i,P)),f)
# end

# function =={P,M}(f1::FiniteField{P,M}, f2::FiniteField{P,M})
#   return f1.polynomialModulus == f2.polynomialModulus
# end

# function show(io::IO, e::FiniteFieldElement)
#   print(io, e.element, " \u2208 ", e.f)
# end

# function =={P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   return e1.element == e2.element && e1.f == e2.f
# end

# function +{P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   @assert e1.f == e2.f
#   return FiniteFieldElement{P,M}(e1.element+e2.element,e1.f)
# end

# function -{P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   @assert e1.f == e2.f
#   return FiniteFieldElement{P,M}(e1.element-e2.element,e1.f)
# end

# function *{P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   @assert e1.f == e2.f
#   return FiniteFieldElement{P,M}(e1.element*e2.element,e1.f)
# end

# function divrem{P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   @assert e1.f == e2.f
#   q,r = divrem(e1.element,e2.element)
#   return (FiniteFieldElement{P,M}(q,e1.f), FiniteFieldElement{P,M}(r,e1.f))
# end

# function div{P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   q,_ = divrem(e1,e2)
#   return q
# end

# function rem{P,M}(e1::FiniteFieldElement{P,M}, e2::FiniteFieldElement{P,M})
#   q,r = divrem(e1,e2)
#   return r
# end

# function -{P,M}(e::FiniteFieldElement{P,M})
#   return FiniteFieldElement{P,M}(-e.element,e1.f)
# end

# function abs{P,M}(e::FiniteFieldElement{P,M})
#   return abs(e.element)
# end

# function ^{P,M}(e::FiniteFieldElement{P,M}, n::Integer)
#   return FiniteFieldElement{P,M}(e.element^n, e.f)
# end

# function inverse{P,M}(e::FiniteFieldElement{P,M})
#   if e == FiniteFieldElement{P,M}(Polynomial([makeModular(0,P)]),e.f)
#     error("Divide by zero")
#   end

#   x,y,d = extendedEuclideanAlgorithm(e.element, e.f.polynomialModulus)
#   return FiniteFieldElement{P,M}(x, e.f) * FiniteFieldElement{P,M}(Polynomial([inverse(d.coefficients[1])]), e.f)
# end


text:
@testset "$TEST $G" begin

g3 = PathGraph(5, G)
d = EdgeMap(g3, [ 0 1 2 3 4; 5 0 6 7 8; 9 10 0 11 12; 13 14 15 0 16; 17 18 19 20 0])
z = floyd_warshall_shortest_paths(g3, d)
@test z.dists[3,:][:] == [7, 6, 0, 11, 27]
@test z.parents[3,:][:] == [2, 3, 0, 3, 4]

@test enumerate_paths(z)[2][2] == []
@test enumerate_paths(z)[2][4] == enumerate_paths(z,2)[4] == enumerate_paths(z,2,4) == [2,3,4]

end # testset


text:
# Parser (SLR(1)) template

parser_pos_start = pos_new()
parser_pos_end = pos_new()

tok2index = Dict()
for (i, t) in enumerate(tokens)
    tok2index[t] = i
end

# Returns (start, end) describing the start / end
# location of the current token being produced
function parser_pos()
    global parser_pos_start, parser_pos_end

    return (parser_pos_start, parser_pos_end)
end

function _pparse(table, tokens, prods, tok2index, produce_rules)
    global parser_pos_start, parser_pos_end

    stack = Array{Any}([1])
    while true
        @assert length(stack) > 0 && length(tokens) > 0 "Unexpected end of file"

        state = stack[end]
        token = tokens[1]
        action = table[state, tok2index[token]]

        if action == nothing
            if token.id == "\$"
                error("Syntax error : Unexpected end of file")
            else
                error("Syntax error from $(token.start_pos) to $(token.end_pos) : Invalid token $token (data = $(token.data))")
            end
        elseif action == "acc"
            @assert length(tokens) == 1 "Unexpected end of file"

            # Reduce using the first rule
            right = stack[2:2:end]
            rule = prods[begin]

            new_tok = deepcopy(rule.left)
            new_tok.start_pos = right[begin].start_pos
            new_tok.end_pos = right[end].end_pos

            parser_pos_start = deepcopy(new_tok.start_pos)
            parser_pos_end = deepcopy(new_tok.end_pos)

            right = [r.data for r in right]
            try
                new_tok.data = produce_rules[begin](right...)
            catch e
                error("Can't produce token $new_tok, invalid rule (got $(length(right)) right tokens)")
            end

            # Return the root of the ast
            return new_tok.data
        elseif action[1] == 'R'
            # Reduce
            i = Base.parse(Int, action[2:end])
            rule = prods[i]

            # Pop right
            right_toks = stack[length(stack) - length(rule.right) * 2 + 1:2:length(stack)]
            splice!(stack, length(stack) - length(rule.right) * 2 + 1:length(stack))

            # State when we start to produce this token
            old_state = stack[end]

            # Push left
            idx = tok2index[rule.left]

            goto = table[old_state, idx]
            @assert goto != nothing && goto[1] == 'S' "Can't go to new state after using rule #$i (producing $(rule.left) with $(rule.right))"

            new_state = Base.parse(Int, goto[2:end])

            new_tok = deepcopy(rule.left)

            new_tok.start_pos = right_toks[begin].start_pos
            new_tok.end_pos = right_toks[end].end_pos
            parser_pos_start = deepcopy(new_tok.start_pos)
            parser_pos_end = deepcopy(new_tok.end_pos)

            right = [r.data for r in right_toks]
            try
                new_tok.data = produce_rules[i](right...)
            catch e
                error("Can't produce token $new_tok, invalid rule (got $(length(right)) right tokens)")
            end

            push!(stack, new_tok)
            push!(stack, new_state)
        elseif action[1] == 'S'
            # Shift
            next_state = Base.parse(Int, action[2:end])

            # Accept the token, push it on the stack and the new state
            popfirst!(tokens)
            push!(stack, token)
            push!(stack, next_state)
        end
    end
end

# Parser parse
function pparse(tokens)
    global tok2index
    global prods
    global table
    global produce_rules

    return _pparse(table, tokens, prods, tok2index, produce_rules)
end


text:
function lorenz()
    dt = 0.01
    n  = 20000000

    a = 10.
    b = 28.
    c = 8/3

    x = zeros(n)
    y = zeros(n)
    z = zeros(n)

    x[1] = 0
    y[1] = 1.
    z[1] = 1.05

    for i = 1:n-1
        x[i+1] = x[i] + dt * (- a    * x[i] + a * y[i])
        y[i+1] = y[i] + dt * (- x[i] * z[i] + b * x[i] - y[i])
        z[i+1] = z[i] + dt * ( x[i]  * y[i] - c * z[i] )
    end
    return x,y,z
end

gx,gy,gz = lorenz()

max_stars_repo_name: Chachay/JuPySwTest

max_issues_repo_name: Chachay/JuPySwTest

max_forks_repo_name: Chachay/JuPySwTest