text:
# Map data types from Julia to DNNL
dnnl_type(::Type{Float16}) = Lib.dnnl_f16
dnnl_type(::Type{Float32}) = Lib.dnnl_f32
dnnl_type(::Type{Int32}) = Lib.dnnl_s32
dnnl_type(::Type{Int8}) = Lib.dnnl_s8
dnnl_type(::Type{UInt8}) = Lib.dnnl_u8

# Try to forward Numbers to their underlying type
dnnl_type(x::T) where {T<:Number} = dnnl_type(T)
dnnl_type(::Type{T}) where {T} = error("No DNNL type for type $T")
dnnl_type(::T) where {T} = error("No DNNL type for $T")

# Automatic argument conversion.
dnnl_convert(::Type{T}) where {T} = dnnl_type(T)

# Make a DIMS array
# NOTE: The OneDNN C-API expects a pointer, so we can't just pass a tuple.
# We either need to pass an array, or a Ref{Tuple}.
# Hwere, we choose to do the latter.
function dnnl_dims(x::NTuple{N,Int64}) where {N}
    # Need to reverse the order of dimensions because C/C++ is row-major while Julia
    # is column major.
    #
    # This leads to all kinds of translation schenanigans ...
    f(i) = (i <= length(x)) ? Lib.dnnl_dim_t(x[i]) : zero(Lib.dnnl_dim_t)
    return ntuple(i -> f(i), Val(Lib.DNNL_MAX_NDIMS))
end
dnnl_dims(x::NTuple{Lib.DNNL_MAX_NDIMS,Int64}) = x

dnnl_dims(x::AbstractArray) = dnnl_dims(strides(x))
dnnl_dims() = ntuple(i -> zero(Int64), Val(Lib.DNNL_MAX_NDIMS))
dnnl_dims(::Tuple{}) = dnnl_dims()

# Automatically wrap in a "Ref" when calling an API function.
# This is almost always passed by reference.
# If it's ever passed by value, then we'll have to add an escape hatch.
dnnl_convert(x::NTuple{N,Int64}) where {N} = Ref(dnnl_dims(x))

# Formats
default_strides(size::Tuple{Vararg{Int,N}}) where {N} = Base.size_to_strides(1, size...)
dnnl_format_any() = Lib.dnnl_format_tag_any

#####
##### Memory Desc
#####

const MemoryDesc = Lib.dnnl_memory_desc_t
dnnl_convert(x::NTuple{N,MemoryDesc}) where {N} = Ref(x)

memorydesc(x::MemoryDesc) = x
function Base.cconvert(::Type{Ptr{Lib.dnnl_memory_desc_t}}, x::MemoryDesc)
    return Base.cconvert(Ptr{Lib.dnnl_memory_desc_t}, Ref(x))
end

# Specialize on typed dimensions to make typestable
logicalsize(md::MemoryDesc) = reverse(md.dims[1:min(12, md.ndims)])
logicalsize(md::MemoryDesc, v::Val) = reverse(ntuple(i -> md.dims[i], v))
function Base.strides(md::MemoryDesc, v::Val)
    return reverse(ntuple(i -> md.format_desc.strides[i], v))
end

function padded_size(md::MemoryDesc, v::Val)
    return reverse(ntuple(i -> md.padded_dims[i], v))
end
Base.ndims(md::MemoryDesc) = md.ndims

function Base.show(io::IO, md::MemoryDesc)
    @nospecialize

    # Perform the "min" check to handle corner cases where the memory layout is too exotic
    # for this simple showing method.
    ndims = min(md.ndims, 12)
    size = logicalsize(md)
    data_type = md.data_type
    padded_dims = reverse(md.padded_dims[1:ndims])
    format_kind = md.format_kind
    padded_offsets = md.padded_offsets[1:ndims]

    # Additional information if the format is "blocked"
    extra_format = ""
    if format_kind == Lib.dnnl_blocked
        #blocking_desc = md.format_desc.blocking
        blocking_desc = md.format_desc
        num_inner_blocks = min(12, blocking_desc.inner_nblks)
        extra_format = """
                strides: $(blocking_desc.strides[1:ndims])
                num inner blocks: $(num_inner_blocks)
                inner blocks: $(reverse(blocking_desc.inner_blks[1:num_inner_blocks]))
                inner indexes: $(reverse(blocking_desc.inner_idxs[1:num_inner_blocks]))
        """
    end

    str = """
    OneDNN Memory Description
        ndims: $ndims
        size: $size
        datatype: $data_type
        format kind: $format_kind
    $extra_format
        padded dims: $padded_dims
        padded offsets: $padded_offsets
    """
    print(io, str)
    return nothing
end

#####
##### Memory Desc constructors in all their glory!
#####

memorydesc(x::AbstractArray{T}) where {T} = memorydesc(T, size(x), strides(x))
function memorydesc(
    datatype, dims::NTuple{N,Int}, strides::NTuple{N,Int} = default_strides(dims)
) where {N}
    handle = Ref{MemoryDesc}()
    @apicall dnnl_memory_desc_init_by_strides(
        handle, N, reverse(dims), datatype, reverse(strides)
    )
    return unwrap_ref(handle)
end

# convenience creation by tag.
function memorydesc(
    datatype, dims::NTuple{N,Int}, tag::Union{Lib.dnnl_format_tag_t,UInt32}
) where {T,N}
    handle = Ref{MemoryDesc}()
    @apicall dnnl_memory_desc_init_by_tag(handle, N, reverse(dims), datatype, tag)
    return unwrap_ref(handle)
end

toany(a::MemoryDesc) = memorydesc(a.data_type, logicalsize(a), dnnl_format_any())

isany(a::Ptr{MemoryDesc}) = isany(unsafe_load(a))
isany(a::MemoryDesc) = a.format_kind == Lib.dnnl_format_kind_any

function Base.:(==)(a::MaybeRef{MemoryDesc}, b::MaybeRef{MemoryDesc})
    return Bool(Lib.dnnl_memory_desc_equal(wrap_ref(a), wrap_ref(b)))
end

getbytes(a::MaybeRef{MemoryDesc}) = signed(Lib.dnnl_memory_desc_get_size(wrap_ref(a)))

#####
##### Memory
#####

# Often, we don't want to specialize on the layout of the array, instead allowing oneDNN
# to work on its own
struct Memory{T,N,A<:AbstractArray{T}} <: AbstractArray{T,N}
    # The underlying array that is supplying the data.
    array::A
    offset::Int

    # Keep around some information about size and padding.
    logicalsize::NTuple{N,Int}

    # Memory object from DNNL
    memory::MemoryPtr
end

function Memory(
    array::A, offset, dims::NTuple{N,Int}, memory::MemoryPtr
) where {T,N,A<:AbstractArray{T}}
    return Memory{T,N,A}(array, convert(Int64, offset), dims, memory)
end

function Base.convert(::Type{Memory{T,N,A}}, x::Memory{T,N,B}) where {T,N,A,B}
    return Memory(convert(A, x.array), x.offset, x.logicalsize, x.memory)
end

Base.sizeof(x::Memory) = getbytes(memorydesc_ptr(x))
toany(x::Memory) = toany(memorydesc(x))

logicalsize(x::Memory) = size(x)
Base.strides(x::Memory{T,N}) where {T,N} = strides(memorydesc(x), Val(N))
padded_size(x::Memory{T,N}) where {T,N} = padded_size(memorydesc(x), Val(N))

Base.parent(x::Memory) = x.array
function ChainRulesCore.rrule(::typeof(Base.parent), x::Memory)
    return parent(x), Δ -> (ChainRulesCore.NoTangent(), Δ)
end

arraytype(::Memory{T,N,A}) where {T,N,A} = A

function Base.show(io::IO, x::Memory)
    print(io, "Opaque Memory $(logicalsize(x))")
    x.offset != 1 && print(io, " - SubArray")
end
Base.show(io::IO, ::MIME"text/plain", x::Memory) = show(io, x)

# for creating OneDNN arguments
@inline access_pointer(x, offset, context) = pointer(x, offset)
function setptr!(x::Memory{T}, context::AccessContext = Reading()) where {T}
    ptr = access_pointer(x.array, x.offset, context)
    @apicall dnnl_memory_set_data_handle_v2(x.memory, ptr, global_stream())
end

function Base.cconvert(
    ::Type{T}, x::Memory
) where {T<:Union{Lib.dnnl_memory_t,Ptr{Lib.dnnl_memory_t}}}
    setptr!(x)
    return Base.cconvert(T, x.memory)
end

Base.cconvert(::Type{Ptr{Lib.dnnl_memory_desc_t}}, x::Memory) = memorydesc_ptr(x)

Base.elsize(::Type{<:Memory{T}}) where {T} = sizeof(T)
function Base.unsafe_convert(::Type{Ptr{T}}, x::Memory{T}) where {T}
    return pointer(x.array)
end

# For constructing DNNL arguments.
function dnnl_exec_arg(x::Memory, context::AccessContext = Reading())
    setptr!(x, context)
    return x.memory
end

# Try to remove as many layers of wrapping around `A` as possible.
# Since all of the dimension and layout information will be stored in the OneDNN
# `memorydesc`, we don't need to hold onto it on the Julia level, which can potentially
# cause down-stream type instabilities.
function Memory(A::AbstractArray)
    return Memory(ancestor(A), _offset(A), size(A), MemoryPtr(A))
end

_offset(::AbstractArray) = one(Int64)
_offset(x::Base.SubArray) = Base.first_index(x)

Memory(M::Memory) = M

function ChainRulesCore.rrule(::Type{<:Memory}, x)
    return (Memory(x), Δ -> (ChainRulesCore.NoTangent(), Δ))
end

# Convenience method for creating destination memories from a source memory.
Base.size(M::Memory) = M.logicalsize
Base.eltype(M::Memory{T}) where {T} = T

function Base.getindex(::Memory, I::Vararg{Int,N}) where {N}
    return error("Cannot index opaque memory formats")
end

function Base.setindex!(::Memory, v, I::Vararg{Int,N}) where {N}
    return error("Cannot index opaque memory formats")
end

memorydesc(M::Memory) = unsafe_load(memorydesc_ptr(M))
function memorydesc_ptr(M::Memory)
    md = Ref{Ptr{MemoryDesc}}()
    @apicall Lib.dnnl_memory_get_memory_desc(M.memory, md)
    return unwrap_ref(md)
end

#####
##### Lazy Transpose
#####

# General idea: swap the dims and strides.
# TODO: Need to validate that this is a blocked layout with no tiling ...
function Base.adjoint(M::Memory{T,2}) where {T}
    dims = size(M)
    strides = Base.strides(memorydesc(M), Val(2))

    reversed_dims = reverse(dims)
    desc = memorydesc(T, reversed_dims, reverse(strides))
    memory = MemoryPtr(parent(M), desc)
    return Memory(parent(M), M.offset, reversed_dims, memory)
end

function Base.permutedims(M::Memory{T,N}, perm::NTuple{N,Int}) where {T,N}
    dims = size(M)
    strides = Base.strides(memorydesc(M), Val(N))
    dims_permuted = unsafe_permute(dims, perm)
    strides_permuted = unsafe_permute(strides, perm)

    desc = memorydesc(T, dims_permuted, strides_permuted)
    memory = MemoryPtr(parent(M), desc)
    return Memory(parent(M), M.offset, dims_permuted, memory)
end

function unsafe_permute(a::NTuple{N,Int}, b::NTuple{N,Int}) where {N}
    return ntuple(i -> @inbounds(a[@inbounds b[i]]), Val(N))
end

#####
##### Construct more memories!!
#####

function Base.similar(
    x::Memory{U,M},
    ::Type{T} = eltype(x),
    dims::NTuple{N,Int} = size(x),
    desc::MemoryDesc = (M == N && U === T) ? memorydesc(x) : memorydesc(T, dims),
) where {U,T,M,N}
    # Number of bytes to allocate.
    # Since OneDNN is free to reorder and pad, we need to explicitly ask it.
    bytes = getbytes(desc)

    # Allocate the output array.
    # This will be allocated as just a plain vector with dimensions padded with ones so it
    # has the same dimension as the wrapped "Memory"
    padded_dims = (div(bytes, sizeof(T)), ntuple(_ -> 1, Val(N - 1))...)
    out = similar(x.array, T, padded_dims)

    # Since we specifically created this array, the offset will always start atone.
    return Memory(out, 1, dims, MemoryPtr(out, desc))
end

materialize(x::AbstractArray, args...; kw...) = x
function materialize(M::Memory{T,N}; allowreorder = true) where {T,N}
    # Check if this memory is already in the requested layout.
    # If so, return the underlying array.
    desired_strides = default_strides(logicalsize(M))
    actual_strides = strides(M)

    # In order to return the underlying object, we need to ensure that:
    # 1. The length of the wrapped object is the same as the length of the Memory.
    # This helps handle views correctly.
    #
    # 2. Strides are the same[
    if length(parent(M)) == length(M) && desired_strides == actual_strides
        return reshape(parent(M), logicalsize(M))
    end

    if !allowreorder
        msg = """
        Expected strides: $desired_strides.
        Found strides: $actual_strides.
        """
        throw(ArgumentError(msg))
    end

    desc = memorydesc(T, logicalsize(M), desired_strides)
    return reshape(parent(reorder(desc, M)), logicalsize(M))
end

function ChainRulesCore.rrule(
    ::typeof(materialize), x, args::Vararg{Any,N}; kw...
) where {N}
    return materialize(x, args...; kw...),
    Δ -> (ChainRulesCore.NoTangent(), Δ, ntuple(_ -> ChainRulesCore.NoTangent(), Val(N)))
end

#####
##### Reshape
#####

function Base.reshape(memory::Memory{T}, dims::NTuple{N,Int}) where {T,N}
    md = Ref{MemoryDesc}()
    @apicall dnnl_memory_desc_reshape(md, memory, N, Ref(reverse(dims)))
    new_memory = MemoryPtr(parent(memory), md)
    return Memory(parent(memory), memory.offset, dims, new_memory)
end

#####
##### Bridge to TiledArrays
#####

function generate_linear_indices(memory::Memory{T,N}) where {T,N}
    md = OneDNN.memorydesc(memory)
    _layout = layout(md)
    _size = logicalsize(md, Val(N))
    _padded_size = padded_size(md, Val(N))
    return generate_linear_indices(_layout, _size, _padded_size)
end


text:
for (depfun, newfun) in Dict(
    :num_columns => :ncol,
    :num_rows => :nrow,
    :bind_value! => :bind!,
    :bind_variable! => :bind!,
    :execute! => :execute,
    :close! => :close,
    :commit! => :commit,
    :rollback! => :rollback,
    :fetch! => :fetch,
    :fetch_row! => :fetchrow,
    :fetch_rows! => :fetchrows

)
    @eval begin
        ($depfun)(args...) = error($depfun, " was renamed to ", $newfun, ". Check release notes at https://github.com/felipenoris/Oracle.jl/releases.")
    end
end

ext: jl


text:
unitrng(; eps = 1e-5) = 0 + eps:0.005:1 - eps

function fx(xrv::RandVar;
            ΩT = defΩ(),
            xdim = 1,
            verbose = false)
  ω = ΩT()
  xrv(ω)
  function f(x)
    ω = Space.update(ω, xdim, x)
    y = xrv(ω)
    verbose && println(x, " -> ", y)
    y
  end
end

function fxy(xrv::RandVar;
                  ΩT = defΩ(),
                  xdim = 1,
                  ydim = 2)
  ω = ΩT()
  xrv(ω)
  function f(x, y)
    ω = Space.update(ω, xdim, x)
    ω = Space.update(ω, ydim, y)
    xrv(ω)
  end
end

"""
Contour Plot of two dimensions of Ω

```julia
using Omega
x = normal(0.0, 1.0)
y = normal(0.0, 1.0)
c1 = logerr(x ==ₛ y)
c2 = logerr(x >ₛ y)
ωcontour(c2)
```

"""
function ωcontour(xrv::RandVar;
                  ΩT = defΩ(),
                  xdim = 1,
                  ydim = 2,
                  xrng = unitrng(),
                  yrng = unitrng(),
                  plt = Plots.plot(),
                  fill = true,
                  aspect_ratio = 1,
                  kwargs...)
  ω = ΩT()
  xrv(ω)
  function f(x, y)
    ω = Space.update(ω, xdim, x)
    ω = Space.update(ω, ydim, y)
    @show xrv(ω)
  end
  Plots.contour!(plt, xrng, yrng, f; fill = fill, aspect_ratio = aspect_ratio, kwargs...)
end

function ωcontourmk(xrv::RandVar;
                  ΩT = defΩ(),
                  xdim = 1,
                  ydim = 2,
                  xrng = unitrng(),
                  yrng = unitrng(),
                  aspect_ratio = 1,
                  kwargs...)
  ω = ΩT()
  xrv(ω)
  function f(x, y)
    ω = Space.update(ω, xdim, x)
    ω = Space.update(ω, ydim, y)
    @show abs(xrv(ω)) / 20000
  end
  AbstractPlotting.surface(xrng, yrng, f; kwargs...)
end

# "x ∈ [0, 1]?"
# isunit(x) = 0.0 <= x <= 1.0

"Plot path in Omega Space from sequence of samples"
function ωtrace(ωs::Vector;
                   xdim = 1,
                   ydim = 2,
                   plt = plot())
  xs, ys = ntranspose(map(ω -> (getdim(ω, xdim), getdim(ω, ydim)), ωs))
  plot!(plt, xs, ys, arrow = :arrow, linealpha = 0.5, legend = false)
end


text:
nx_add = 1400000
nx_new = 36000
pref = path * "/SampleFiles/"
cfile = pref * "Restricted/03_02_27_20140927.sjis.ch"
files = String[ "UW/00012502123W"               "uw"            "_"
                "SEGY/03.334.12.09.00.0362.1"   "passcal"       "pa-full"
                "SEGY/03.334.12.09.00.0362.1"   "passcal"       "passcal"
                "SAC/test_le.sac"               "sac"           "_"
                "SEED/test.mseed"               "mseed"         "_"
                "Restricted/2014092709*.cnt"    "win32"         "win"
                "ASCII/0215162000.c00"          "lennartz"      "_"
                "ASCII/geo-slist.csv"           "geocsv.slist"  "_"
                "Restricted/SHW.UW.mseed"       "mseed"         "lo-mem"
                "Restricted/test_rev_1.segy"    "segy"          "full"
                "Restricted/test_rev_1.segy"    "segy"          "_"
                "SAC/test_be.sac"               "sac"           "full"
                "SAC/test_be.sac"               "sac"           "_"
                "ASCII/geo-tspair.csv"          "geocsv"        "_"
                ]

# Do not run on Appveyor; it can't access restricted files, so this breaks
printstyled("  read_data\n", color=:light_green)
nf = size(files,1)
for n = 1:nf
  fname = pref * files[n,1]
  if (occursin("Restricted", fname)==false || (has_restricted==true))
    fwild = fname[1:end-1] * "*"
    f_call = files[n,2]
    opt = files[n,3]
    printstyled(string("    ", n, "/", nf, " ", f_call, "\n"), color=:light_green)
    if opt == "pa-full"
      S = verified_read_data("passcal", fname, full=true)
      S = verified_read_data("passcal", fwild, full=true)
      S = read_data("passcal", fwild, full=true, memmap=true)
      S = read_data(fname)
      S = read_data(fname, full=true)
    elseif opt == "win"
      S = verified_read_data(f_call, fwild, cf=cfile)
      S = read_data(f_call, fwild, memmap=true, cf=cfile)
      S = read_data(fwild, cf=cfile)
      S = read_data(f_call, [fwild], cf=cfile)
    elseif opt == "slist"
      S = verified_read_data("geocsv.slist", fname)
      S = verified_read_data("geocsv.slist", fwild)
      S = read_data("geocsv.slist", fwild, memmap=true)
      S = read_data(fwild)
    elseif opt == "lo-mem"
      S = verified_read_data(f_call, fname, nx_new=nx_new, nx_add=nx_add)
      S = verified_read_data(f_call, fwild, nx_new=nx_new, nx_add=nx_add)
      S = read_data(f_call, fwild, nx_new=nx_new, nx_add=nx_add, memmap=true)
    elseif opt == "full"
      S = verified_read_data(f_call, fname, full=true)
      S = verified_read_data(f_call, fwild, full=true)
      S = read_data(f_call, fwild, full=true, memmap=true)
      S = read_data(fwild, full=true)
    else
      S = verified_read_data(f_call, fname)
      S = read_data(f_call, fname, memmap=true)
      if f_call == "uw"
        fwild = fname[1:end-3]*"*"*"W"
      end
      S = read_data(f_call, fwild, memmap=true)
      S = read_data(fwild)
    end
  end
end

# Test for reading a String array of file names
pref = path * "/SampleFiles/"
files = pref .* ["SAC/test_be.sac", "SAC/test_le.sac"]
S = verified_read_data("sac", files, vl=true)
for f in files
  for i in 1:S.n
    @test any([occursin(abspath(f), n) for n in S.notes[i]])
  end
end

# Check that other array methods work
S1 = SeisData()
S2 = SeisData()
read_data!(S1, "sac", files, vl=true)
read_data!(S2, files, vl=true)
S3 = read_data(files, vl=true)
@test S == S1 == S2 == S3

# Does a string array read the same way as a wildcard read?
compare_SeisData(S, verified_read_data("sac", pref .* "SAC/test*sac", vl=true))

# source logging
printstyled("    logging\n", color=:light_green)
uwf1 = joinpath(path, "SampleFiles/UW/99011116541")
uwf4 = joinpath(path, "SampleFiles/UW/00012502123W")
p1 = abspath(uwf1*"W")
p2 = abspath(uwf4)

S = read_data("uw", [uwf1*"W", uwf4])
for i in 1:S.n
  if length(S.t[i]) == 6
    @test S.src[i] == p2
  else
    @test S.src[i] == (S.t[i][1,2] ≥ 946684800000000 ? p2 : p1)
  end
end

@test_throws ErrorException verified_read_data("deez", "nutz.sac")
nothing


text:
mutable struct TempFileWrapper
    file::Union{String, Nothing}
end

# structure for iterating over a csv file
# no automatic type inference is done, but types are allowed to be passed
# for as many columns as desired; `CSV.detect(row, i)` can also be used to
# use the same inference logic used in `CSV.File` for determing a cell's typed value
struct Rows{IO, customtypes, V, stringtype}
    name::String
    names::Vector{Symbol} # only includes "select"ed columns
    columns::Vector{Column}
    columnmap::Vector{Int} # maps "select"ed column index to actual file column index
    buf::IO
    datapos::Int
    datarow::Int
    len::Int
    limit::Int
    options::Parsers.Options
    reusebuffer::Bool
    values::Vector{V} # once values are parsed, put in values; allocated on each iteration if reusebuffer=false
    lookup::Dict{Symbol, Int}
    numwarnings::Base.RefValue{Int}
    maxwarnings::Int
    ctx::Context
    tempfile::TempFileWrapper
end

function Base.show(io::IO, r::Rows)
    println(io, "CSV.Rows(\"$(r.name)\"):")
    println(io, "Size: $(length(r.columns))")
    show(io, Tables.schema(r))
end

"""
    CSV.Rows(source; kwargs...) => CSV.Rows

Read a csv input returning a `CSV.Rows` object.

The [`input`](@ref input) argument can be one of:
  * filename given as a string or FilePaths.jl type
  * a `Vector{UInt8}` or `SubArray{UInt8, 1, Vector{UInt8}}` byte buffer
  * a `CodeUnits` object, which wraps a `String`, like `codeunits(str)`
  * a csv-formatted string can also be passed like `IOBuffer(str)`
  * a `Cmd` or other `IO`
  * a gzipped file (or gzipped data in any of the above), which will automatically be decompressed for parsing

To read a csv file from a url, use the HTTP.jl package, where the `HTTP.Response` body can be passed like:
```julia
f = CSV.Rows(HTTP.get(url).body)
```

For other `IO` or `Cmd` inputs, you can pass them like: `f = CSV.Rows(read(obj))`.

While similar to [`CSV.File`](@ref), `CSV.Rows` provides a slightly different interface, the tradeoffs including:
  * Very minimal memory footprint; while iterating, only the current row values are buffered
  * Only provides row access via iteration; to access columns, one can stream the rows into a table type
  * Performs no type inference; each column/cell is essentially treated as `Union{String, Missing}`, users can utilize the performant `Parsers.parse(T, str)` to convert values to a more specific type if needed, or pass types upon construction using the `type` or `types` keyword arguments

Opens the file and uses passed arguments to detect the number of columns, ***but not*** column types (column types default to `String` unless otherwise manually provided).
The returned `CSV.Rows` object supports the [Tables.jl](https://github.com/JuliaData/Tables.jl) interface
and can iterate rows. Each row object supports `propertynames`, `getproperty`, and `getindex` to access individual row values.
Note that duplicate column names will be detected and adjusted to ensure uniqueness (duplicate column name `a` will become `a_1`).
For example, one could iterate over a csv file with column names `a`, `b`, and `c` by doing:

```julia
for row in CSV.Rows(file)
    println("a=\$(row.a), b=\$(row.b), c=\$(row.c)")
end
```

$KEYWORD_DOCS
"""
function Rows(source::ValidSources;
    # file options
    # header can be a row number, range of rows, or actual string vector
    header::Union{Integer, Vector{Symbol}, Vector{String}, AbstractVector{<:Integer}}=1,
    normalizenames::Bool=false,
    # by default, data starts immediately after header or start of file
    datarow::Integer=-1,
    skipto::Integer=-1,
    footerskip::Integer=0,
    transpose::Bool=false,
    comment::Union{String, Nothing}=nothing,
    ignoreemptyrows::Bool=true,
    ignoreemptylines=nothing,
    select=nothing,
    drop=nothing,
    limit::Union{Integer, Nothing}=nothing,
    buffer_in_memory::Bool=false,
    # parsing options
    missingstrings=String[],
    missingstring="",
    delim::Union{Nothing, Char, String}=nothing,
    ignorerepeated::Bool=false,
    quoted::Bool=true,
    quotechar::Union{UInt8, Char}='"',
    openquotechar::Union{UInt8, Char, Nothing}=nothing,
    closequotechar::Union{UInt8, Char, Nothing}=nothing,
    escapechar::Union{UInt8, Char}='"',
    dateformat::Union{String, Dates.DateFormat, Nothing, AbstractDict}=nothing,
    dateformats=nothing,
    decimal::Union{UInt8, Char}=UInt8('.'),
    truestrings::Union{Vector{String}, Nothing}=TRUE_STRINGS,
    falsestrings::Union{Vector{String}, Nothing}=FALSE_STRINGS,
    # type options
    type=nothing,
    types=nothing,
    typemap::Dict=Dict{Type, Type}(),
    pool=false,
    downcast::Bool=false,
    stringtype::StringTypes=PosLenString,
    lazystrings::Bool=stringtype === PosLenString,
    strict::Bool=false,
    silencewarnings::Bool=false,
    maxwarnings::Int=100,
    debug::Bool=false,
    parsingdebug::Bool=false,
    validate::Bool=true,
    reusebuffer::Bool=false,
    )
    ctx = @refargs Context(source, header, normalizenames, datarow, skipto, footerskip, transpose, comment, ignoreemptyrows, ignoreemptylines, select, drop, limit, buffer_in_memory, nothing, nothing, nothing, 0, nothing, missingstrings, missingstring, delim, ignorerepeated, quoted, quotechar, openquotechar, closequotechar, escapechar, dateformat, dateformats, decimal, truestrings, falsestrings, type, types, typemap, pool, downcast, lazystrings, stringtype, strict, silencewarnings, maxwarnings, debug, parsingdebug, validate, true)
    foreach(col -> col.pool = 0.0, ctx.columns)
    allocate!(ctx.columns, 1)
    values = all(x->x.type === ctx.stringtype && x.anymissing, ctx.columns) && ctx.stringtype === PosLenString ? Vector{PosLen}(undef, ctx.cols) : Vector{Any}(undef, ctx.cols)
    columnmap = collect(1:ctx.cols)
    for i = ctx.cols:-1:1
        col = ctx.columns[i]
        if col.willdrop
            deleteat!(ctx.names, i)
            deleteat!(columnmap, i)
        end
    end
    lookup = Dict(nm=>i for (i, nm) in enumerate(ctx.names))
    tempfile = TempFileWrapper(ctx.tempfile)
    if tempfile.file !== nothing
        finalizer(tempfile) do x
            rm(x.file; force=true)
        end
    end
    return Rows{typeof(ctx.buf), ctx.customtypes, eltype(values), ctx.stringtype}(
        ctx.name,
        ctx.names,
        ctx.columns,
        columnmap,
        ctx.buf,
        ctx.datapos,
        ctx.datarow,
        ctx.len,
        ctx.limit,
        ctx.options,
        reusebuffer,
        values,
        lookup,
        Ref(0),
        ctx.maxwarnings,
        ctx,
        tempfile
    )
end

Tables.isrowtable(::Type{<:Rows}) = true
Tables.schema(r::Rows) = Tables.Schema(r.names, [coltype(x) for x in view(r.columns, r.columnmap)])
Base.eltype(::Rows) = Row2
Base.IteratorSize(::Type{<:Rows}) = Base.SizeUnknown()

@inline function setcustom!(::Type{customtypes}, values, columns, i, setmissing) where {customtypes}
    if @generated
        block = Expr(:block)
        push!(block.args, quote
            error("CSV.jl code-generation error, unexpected column type: $(typeof(column))")
        end)
        for i = 1:fieldcount(customtypes)
            T = fieldtype(customtypes, i)
            vT = vectype(T)
            pushfirst!(block.args, quote
                column = columns[i].column
                if column isa $vT
                    if setmissing
                        setmissing!(column, 1)
                    else
                        @inbounds values[i] = column[1]
                    end
                    return
                end
            end)
        end
        pushfirst!(block.args, Expr(:meta, :inline))
        # @show block
        return block
    else
        # println("generated function failed")
        @inbounds values[i] = columns[i].column[1]
        return
    end
end

function checkwidencolumns!(r::Rows{ct, V}, cols) where {ct, V}
    if cols > length(r.values)
        # we widened while parsing this row, need to widen other supporting objects
        for i = (length(r.values) + 1):cols
            push!(r.values, V === Any ? missing : Base.bitcast(PosLen, Parsers.MISSING_BIT))
            nm = Symbol(:Column, i)
            push!(r.names, nm)
            r.lookup[nm] = length(r.values)
            push!(r.columnmap, i)
        end
    end
    return
end

macro unrollcolumns(setmissing, values, ex)
    return esc(quote
        if column isa MissingVector
            if !($setmissing)
                @inbounds $(values)[i] = missing
            end
        elseif column isa Vector{PosLen}
            $ex
        elseif column isa Vector{Union{Missing, Int8}}
            $ex
        elseif column isa Vector{Union{Missing, Int16}}
            $ex
        elseif column isa Vector{Union{Missing, Int32}}
            $ex
        elseif column isa SVec{Int32}
            $ex
        elseif column isa SVec{Int64}
            $ex
        elseif column isa SVec{Int128}
            $ex
        elseif column isa SVec{Float16}
            $ex
        elseif column isa SVec{Float32}
            $ex
        elseif column isa SVec{Float64}
            $ex
        elseif column isa SVec{InlineString1}
            $ex
        elseif column isa SVec{InlineString3}
            $ex
        elseif column isa SVec{InlineString7}
            $ex
        elseif column isa SVec{InlineString15}
            $ex
        elseif column isa SVec{InlineString31}
            $ex
        elseif column isa SVec{InlineString63}
            $ex
        elseif column isa SVec{InlineString127}
            $ex
        elseif column isa SVec{InlineString255}
            $ex
        elseif column isa SVec2{String}
            $ex
        elseif column isa SVec{Date}
            $ex
        elseif column isa SVec{DateTime}
            $ex
        elseif column isa SVec{Time}
            $ex
        elseif column isa Vector{Union{Missing, Bool}}
            $ex
        elseif column isa Vector{UInt32}
            $ex
        elseif customtypes !== Tuple{}
            setcustom!(customtypes, $values, columns, i, $setmissing)
        else
            error("bad array type: $(typeof(column))")
        end
    end)
end

@inline function Base.iterate(r::Rows{IO, customtypes, V, stringtype}, (pos, len, row)=(r.datapos, r.len, 1)) where {IO, customtypes, V, stringtype}
    (pos > len || row > r.limit) && return nothing
    columns = r.columns
    cols = length(columns)
    for i = 1:cols
        @inbounds column = columns[i].column
        @unrollcolumns true nothing begin
            setmissing!(column, 1)
        end
    end
    pos = parserow(1, 1, r.numwarnings, r.ctx, r.buf, pos, len, 1, r.datarow + row - 2, columns, customtypes)
    cols = length(columns)
    values = r.reusebuffer ? r.values : Vector{V}(undef, cols)
    checkwidencolumns!(r, cols)
    for i = 1:cols
        @inbounds column = columns[i].column
        @unrollcolumns false values begin
            @inbounds values[i] = column[1]
        end
    end
    return Row2{V, stringtype}(r.names, r.columns, r.columnmap, r.lookup, values, r.buf), (pos, len, row + 1)
end

struct Row2{V, stringtype} <: Tables.AbstractRow
    names::Vector{Symbol}
    columns::Vector{Column}
    columnmap::Vector{Int}
    lookup::Dict{Symbol, Int}
    values::Vector{V}
    buf::Vector{UInt8}
end

getnames(r::Row2) = getfield(r, :names)
getcolumns(r::Row2) = getfield(r, :columns)
getcolumnmap(r::Row2) = getfield(r, :columnmap)
getlookup(r::Row2) = getfield(r, :lookup)
getvalues(r::Row2) = getfield(r, :values)
getbuf(r::Row2) = getfield(r, :buf)
getV(::Row2{V}) where {V} = V
getstringtype(::Row2{V, stringtype}) where {V, stringtype} = stringtype

Tables.columnnames(r::Row2) = getnames(r)

Base.checkbounds(r::Row2, i) = 0 < i < length(r)

Tables.getcolumn(r::Row2, nm::Symbol) = Tables.getcolumn(r, getlookup(r)[nm])
Tables.getcolumn(r::Row2, i::Int) = Tables.getcolumn(r, coltype(getcolumns(r)[i]), i, getnames(r)[i])

Base.@propagate_inbounds function Tables.getcolumn(r::Row2, ::Type{T}, i::Int, nm::Symbol) where {T}
    @boundscheck checkbounds(r, i)
    j = getcolumnmap(r)[i]
    values = getvalues(r)
    V = getV(r)
    @inbounds val = j > length(values) ? (V === PosLen ? Parsers.MISSING_BIT : missing) : values[j]
    stringtype = getstringtype(r)
    if V === PosLen
        # column type must be stringtype
        # @show T, stringtype
        @assert T === Union{stringtype, Missing}
        e = getcolumns(r)[j].options.e
        if (val isa PosLen && val.missingvalue) || val == Parsers.MISSING_BIT
            return missing
        elseif stringtype === PosLenString
            return PosLenString(getbuf(r), val, e)
        elseif stringtype === String
            return Parsers.getstring(getbuf(r), val, e)
        end
    else
        # at least some column types were manually provided
        if val isa PosLen
            if val.missingvalue
                return missing
            else
                e = getcolumns(r)[j].options.e
                return PosLenString(getbuf(r), val, e)
            end
        else
            return val
        end
    end
end

@noinline stringsonly() = error("Parsers.parse only allowed on String column types")

Base.@propagate_inbounds function Parsers.parse(::Type{T}, r::Row2, i::Int) where {T}
    @boundscheck checkbounds(r, i)
    @inbounds begin
        j = getcolumnmap(r)[i]
        col = getcolumns(r)[j]
        col.type isa StringTypes || stringsonly()
        poslen = getvalues(r)[j]
        poslen.missingvalue && return missing
        pos = poslen.pos
        res = Parsers.xparse(T, getbuf(r), pos, pos + poslen.len, col.options)
    end
    return Parsers.ok(res.code) ? (res.val::T) : missing
end

Base.@propagate_inbounds function detect(r::Row2, i::Int)
    @boundscheck checkbounds(r, i)
    @inbounds begin
        j = getcolumnmap(r)[i]
        col = getcolumns(r)[j]
        col.type isa StringTypes || stringsonly()
        poslen = getvalues(r)[j]
        poslen.missingvalue && return missing
        pos = poslen.pos
        code, tlen, x, xT = detect(pass, getbuf(r), pos, pos + poslen.len - 1, col.options)
        return x === nothing ? r[i] : x
    end
end

function Parsers.parse(::Type{T}, r::Row2, nm::Symbol) where {T}
    @inbounds x = Parsers.parse(T, r, getlookup(r)[nm])
    return x
end

function detect(r::Row2, nm::Symbol)
    @inbounds x = detect(r, getlookup(r)[nm])
    return x
end


text:
function simurec(
    model::Union{elmod2d,acmod2d},
    sou::Array{MTsource},
    rec::receiver,
    OptCpnt::Array{String},
    RecDataPath::String)

    nT = model.medium.nT
    fid = open(RecDataPath, "a+")
    OptCpnt = checkrepeat(OptCpnt)

    # wIB(fid,rec,model.medium.nT,OptCpnt)
    for it = 1 : nT
        for isn in 1:length(sou)
            addmt!(model.wf, sou[isn], it)
        end
        run!(model)
        wDB(fid,model.wf,rec.BDnloc,model.nwf,OptCpnt)
    end
    close(fid)
end


text:
# # Modeling with JuMP

using PowerSystems
using JuMP


text:
@testset "Initial population reports" begin
    mpSim2 = MPsim()
    filePath = joinpath( "base", "baseConfig" )
    fileName = string( filePath, ".xlsx" )

    XLSX.openxlsx( fileName ) do xf
        MPIO.readInitPop( mpSim2, xf, filePath )
    end  # XLSX.openxlsx( fileName ) do xf

    fName = "report/initpop"
    excelInitPopReport( mpSim2, filename=fName )
    excelInitPopAgeReport( mpSim2, 12, :age, filename=fName, overwrite=false )
    excelInitPopAgeReport( mpSim2, 12, :tenure, filename=fName,
        overwrite=false )
    excelInitPopAgeReport( mpSim2, 12, :timeInNode, filename=fName,
        overwrite=false )
end  # @testset "Initial population reports"


text:
const julia_getModuleAt_request_type = JSONRPC.RequestType("julia/getModuleAt", VersionedTextDocumentPositionParams, String)
const julia_getCurrentBlockRange_request_type = JSONRPC.RequestType("julia/getCurrentBlockRange", VersionedTextDocumentPositionParams, Position)
const julia_getDocAt_request_type = JSONRPC.RequestType("julia/getDocAt", VersionedTextDocumentPositionParams, String)
const julia_getDocFromWord_request_type = JSONRPC.RequestType("julia/getDocFromWord", NamedTuple{(:word,),Tuple{String}}, String)


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Tokei_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Tokei")
JLLWrappers.@generate_main_file("Tokei", UUID("3ac119c9-1236-5556-b556-adc8150b0244"))
end  # module Tokei_jll


text:
using RubiSymbolics: syntax_repl_dict, function_repl_dict, identifier

rubi_utils_file = joinpath(@__DIR__, "../Rubi/IntegrationUtilityFunctions.m")
utils_file = joinpath(@__DIR__, "utils.jl")

orig_file_content = open(rubi_utils_file) do rf
    read(rf, String)
end
out_file_content = replace(orig_file_content,
    "lst[[1,1]][[1]]" => "Part(Part(lst, 1, 1), 1)",
    "Reap[u][[2]]" => "Part(Reap(u), 2)",
    "(lst=MonomialFactor[u,x])[[1]]" => "Part(lst=MonomialFactor(u, x), 1)",
    syntax_repl_dict...,
    )
out_file_content = replace(out_file_content,
    function_repl_dict...,
    )
out_file_content = replace(out_file_content,
    "=#*" => "=#",
    "*#=" => "#=",
    r"(?<function>\w+)::usage\h*=[\h\n]*\"(?<doc>.*?)\";"s => s"\"\g<doc>\"\nfunction \g<function> end",
    "FalseQ(false) = true" => "FalseQ(b::Bool) = !b",
    )
out_file_content = replace(out_file_content,
    r"(\n(?!\n|$))"   => "\n"*" "^4,    # indent
    r"^(?!\n)"      => " "^4,
    )
touch(utils_file)
open(utils_file, "w") do f
    write(f, "@util_collect begin\n")
    write(f, out_file_content)
    write(f, "end\n")
end

# when function(x^2+y...) -> place rule macro in front
# macro will parse & collect rules -> all rules to be applied on rhs after matching
# others: just functions to call

# @rule_ Module @utils Int(...) := ...

# @register_symbolic lhs of := and define rule


text:
#!/usr/bin/env julia

using WriteVTK
using Base.Test
import Compat: UTF8String, readstring

tests = ["multiblock.jl",
         "rectilinear.jl",
         "imagedata.jl",
         "structured.jl",
         "unstructured.jl",
         "pvdCollection.jl"]

# Only toggle to generate new checksums, if new tests are added.
OVERWRITE_CHECKSUMS = false
checksums_file = "checksums.sha1"

checksum_list = readstring(checksums_file)

if OVERWRITE_CHECKSUMS
    csio = open(checksums_file, "w")
end

# Run the test scripts.
for test in tests
    println("TEST (first run): " * test)
    @time outfiles = evalfile(test)::Vector{UTF8String}

    # Check that the generated files match the stored checksums.
    for file in outfiles
        sha = readstring(`sha1sum $file`)
        if OVERWRITE_CHECKSUMS
            write(csio, sha)
        else
            # Returns 0:-1 if string is not found.
            cmp = search(checksum_list, sha)
            @test cmp != 0:-1
        end
    end
    println()
end

OVERWRITE_CHECKSUMS && close(csio)

println("="^60 * "\n")

# Run the tests again, just to measure the time and allocations once all the
# functions have already been compiled.
for test in tests
    println("TEST (second run): " * test)
    @time outfiles = evalfile(test)::Vector{UTF8String}
    println()
end


text:
export enumerate

import Base: enumerate
import Base: show

"""
    enumerate()

Creates an enumerate operator, which converts each value emitted by the source
Observable into a tuple of its order number and the value itself.

The enumerate operator is similar to `scan(Tuple{Int, Int}, (d, c) -> (d, c[2] + 1), (0, 0))` (see [`scan`](@ref)).

# Producing

Stream of type `<: Subscribable{Tuple{Int, L}}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 0 for _ in 1:3 ])
subscribe!(source |> enumerate(), logger())
;

# output

[LogActor] Data: (1, 0)
[LogActor] Data: (2, 0)
[LogActor] Data: (3, 0)
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`scan`](@ref), [`map`](@ref), [`logger`](@ref)
"""
enumerate() = EnumerateOperator()

struct EnumerateOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{Tuple{Int, L}}, operator::EnumerateOperator, source) where L
    return proxy(Tuple{Int, L}, source, EnumerateProxy())
end

operator_right(operator::EnumerateOperator, ::Type{L}) where L = Tuple{Int, L}

struct EnumerateProxy <: ActorProxy end

actor_proxy!(::Type{Tuple{Int, L}}, proxy::EnumerateProxy, actor::A) where { L, A } = EnumerateActor{L, A}(actor, EnumerateActorProps(1))

mutable struct EnumerateActorProps
    current :: Int
end

struct EnumerateActor{L, A} <: Actor{L}
    actor :: A
    props :: EnumerateActorProps
end

getcurrent(actor::EnumerateActor)           = actor.props.current
setcurrent!(actor::EnumerateActor, current) = actor.props.current = current

function on_next!(actor::EnumerateActor{L}, data::L) where L
    current = getcurrent(actor)
    setcurrent!(actor, current + 1)
    next!(actor.actor, (current, data))
end

on_error!(actor::EnumerateActor, err) = error!(actor.actor, err)
on_complete!(actor::EnumerateActor)   = complete!(actor.actor)

Base.show(io::IO, ::EnumerateOperator)         = print(io, "EnumerateOperator()")
Base.show(io::IO, ::EnumerateProxy)            = print(io, "EnumerateProxy()")
Base.show(io::IO, ::EnumerateActor{L}) where L = print(io, "EnumerateActor($L)")


text:
"""
```
alt_inflation()
```

Create an instance of the `AltPolicy` proposed in John Taylor's \"Discretion
versus policy rules in practice\" (1993).
"""
function alt_inflation()
    AltPolicy(:alt_inflation, alt_inflation_eqcond, alt_inflation_solve, color = RGB(1., 0., 0.))
end

function alt_inflation_eqcond(m::AbstractDSGEModel)
    # Get equilibrium condition matrices
    Γ0, Γ1, C, Ψ, Π  = eqcond(m)

    eq   = m.equilibrium_conditions
    endo = m.endogenous_states

    # Zero out old policy rule
    Γ0[eq[:eq_mp], :] .= 0
    Γ1[eq[:eq_mp], :] .= 0
    C[eq[:eq_mp]]     = 0
    Ψ[eq[:eq_mp], :]  .= 0
    Π[eq[:eq_mp], :]  .= 0

    # Put in new policy rule
    Γ0[eq[:eq_mp], endo[:π_t]]   = 1
    Γ0[eq[:eq_mp], endo[:π_star_t]] = -1

    return Γ0, Γ1, C, Ψ, Π
end

function alt_inflation_solve(m::AbstractDSGEModel)
    # Get equilibrium condition matrices
    Γ0, Γ1, C, Ψ, Π = alt_inflation_eqcond(m)

    # Solve model
    TTT_gensys, CCC_gensys, RRR_gensys, eu = gensys(Γ0, Γ1, C, Ψ, Π, 1+1e-6)
    #TTT_gensys, CCC_gensys, RRR_gensys, fmat, fwt, ywt, gev, eu, loose = gensys(Γ0, Γ1, C, Ψ, Π, 1+1e-6)
    if !((eu[1] == 1) & (eu[2] == 1))
        throw(GensysError("Gensys does not give existence"))
    end
    TTT_gensys = real(TTT_gensys)
    RRR_gensys = real(RRR_gensys)
    CCC_gensys = reshape(CCC_gensys, size(CCC_gensys, 1))

    # Augment states
    TTT, RRR, CCC = augment_states(m, TTT_gensys, RRR_gensys, CCC_gensys)

    return TTT, RRR, CCC
end


text:
using IterTools, BlockArrays, LinearAlgebra
using Stheno: MeanFunction, Kernel, CrossKernel, AV, blocks, pairwise, LazyPDMat

"""
    unary_map_tests(f, X::AbstractVector)

Consistency tests intended for use with `MeanFunction`s.
"""
function unary_map_tests(f, X::AbstractVector)
    @test map(f, X) isa AbstractVector
    @test length(map(f, X)) == length(X)
    @test map(f, X) ≈ [f(x) for x in X]
    @test map(f, X) ≈ Stheno._map_fallback(f, X)
end
function unary_map_tests(f, X::BlockData)
    @test map(f, X) isa AbstractBlockVector
    @test length(map(f, X)) == length(X)
    @test map(f, X) ≈ BlockVector([map(f, x) for x in blocks(X)])
end

"""
    binary_map_tests(f, X::AbstractVector, X′::AbstractVector)

Consistency tests intended for use with `CrossKernel`s.
"""
function binary_map_tests(f, X::AbstractVector, X′::AbstractVector)
    @test map(f, X, X′) isa AbstractVector
    @test length(map(f, X, X′)) == length(X)
    @test map(f, X, X′) ≈ [f(x, x′) for (x, x′) in zip(X, X′)]
    @test map(f, X, X′) ≈ Stheno._map_fallback(f, X, X′)
end
function binary_map_tests(f, XB::BlockData, XB′::BlockData)
    @test map(f, XB, XB′) isa AbstractBlockVector
    @test length(map(f, XB, XB′)) == length(XB)
    @test map(f, XB, XB′) ≈
        BlockVector([map(f, X, X′) for (X, X′) in zip(blocks(XB), blocks(XB′))])
end

"""
    binary_map_tests(f, X::AbstractVector)

Consistency tests intended for use with `Kernel`s.
"""
function binary_map_tests(f, X::AbstractVector)
    @test map(f, X) isa AbstractVector
    @test length(map(f, X)) == length(X)
    @test map(f, X) ≈ map(f, X, X)
end
function binary_map_tests(f, X::BlockData)
    @test map(f, X) isa AbstractBlockVector
    @test length(map(f, X)) == length(X)
    @test map(f, X) ≈ map(f, X, X)
end

"""
    pairwise_tests(f, X::ADS, X′::ADS)

Consistency tests intended for use with `CrossKernel`s.
"""
function pairwise_tests(f, X::AbstractVector, X′::AbstractVector)
    N, N′ = length(X), length(X′)
    @test pairwise(f, X, X′) isa AbstractMatrix
    @test size(pairwise(f, X, X′)) == (N, N′)
    @test pairwise(f, X, X′) ≈ reshape([f(x, x′) for (x, x′) in Iterators.product(X, X′)], N, N′)
    @test pairwise(f, X, X′) ≈ Stheno._pairwise_fallback(f, X, X′)
end
function pairwise_tests(f, X::BlockData, X′::BlockData)
    N, N′ = length(X), length(X′)
    @test pairwise(f, X, X′) isa AbstractBlockMatrix
    @test size(pairwise(f, X, X′)) == (N, N′)
    @test pairwise(f, X, X′) ==
        BlockMatrix([pairwise(f, x, x′) for x in blocks(X), x′ in blocks(X′)])
end

"""
    pairwise_tests(f, X::ADS)

Consistency tests intended for use with `Kernel`s.
"""
function pairwise_tests(f, X::AbstractVector; rtol=eps())
    @test pairwise(f, X) isa LazyPDMat{T, <:AbstractMatrix{T}} where T
    @test size(pairwise(f, X)) == (length(X), length(X))
    @test isapprox(pairwise(f, X), pairwise(f, X, X); rtol=rtol)
end
function pairwise_tests(f, X::BlockData; rtol=eps())
    @test pairwise(f, X) isa LazyPDMat{T, <:Symmetric{T, <:AbstractBlockMatrix{T}}} where T
    @test size(pairwise(f, X)) == (length(X), length(X))
    @test pairwise(f, X) ==
        BlockMatrix([pairwise(f, x, x′) for x in blocks(X), x′ in blocks(X)])
end

"""
    mean_function_tests(μ::MeanFunction, X::AbstractVector)

Tests that any mean function `μ` should be able to pass.
"""
function mean_function_tests(μ::MeanFunction, X::AbstractVector)
    __mean_function_tests(μ, X)
    mean_function_tests(μ, BlockData([X]))
    mean_function_tests(μ, BlockData([X, X]))
end
function mean_function_tests(μ::MeanFunction, X::BlockData)
    __mean_function_tests(μ, X)
end
function __mean_function_tests(μ::MeanFunction, X::AbstractVector)

    # Test compulsory interface passes.
    @test !(μ(X[1]) isa Nothing)

    @test hasmethod(eachindex, Tuple{typeof(μ)})

    # Test optional interface.
    unary_map_tests(μ, X)
end

"""
    cross_kernel_tests(k::CrossKernel, X::AbstractVector, X′::AbstractVector)

Tests that any kernel `k` should be able to pass. Requires that `length(X0) == length(X1)`
and `length(X0) ≠ length(X2)`.
"""
function cross_kernel_tests(k::CrossKernel, X0::AV, X1::AV, X2::AV)
    __cross_kernel_tests(k, X0, X1, X2)
    cross_kernel_tests(k, BlockData([X0, X0]), BlockData([X1, X1]), BlockData([X2]))
end
function cross_kernel_tests(k::CrossKernel, X0::BlockData, X1::BlockData, X2::BlockData)
    __cross_kernel_tests(k, X0, X1, X2)
end
function __cross_kernel_tests(k::CrossKernel, X0::AV, X1::AV, X2::AV)
    @assert length(X0) == length(X1)
    @assert length(X0) ≠ length(X2)

    @test hasmethod(eachindex, Tuple{typeof(k), Int})

    binary_map_tests(k, X0, X1)
    pairwise_tests(k, X0, X2)
end

"""
    kernel_tests(k::Kernel, X0::AbstractVector, X1::AbstractVector, X2::AbstractVector)

Tests that any kernel `k` should be able to pass. Requires that `length(X0) == length(X1)`
and `length(X0) ≠ length(X2)`.
"""
function kernel_tests(k::Kernel, X0::AV, X1::AV, X2::AV, rtol::Real=eps())
    __kernel_tests(k, X0, X1, X2, rtol)
    kernel_tests(k, BlockData([X0, X0]), BlockData([X1, X1]), BlockData([X2]), rtol)
end
function kernel_tests(k::Kernel, X0::AV, X1::AV, X2::AV, rtol::Real=eps())
    __kernel_tests(k, X0, X1, X2, rtol)
end
function __kernel_tests(k::Kernel, X0::AV, X1::AV, X2::AV, rtol::Real=eps())
    @assert length(X0) == length(X1)
    @assert length(X0) ≠ length(X2)

    @test hasmethod(eachindex, Tuple{typeof(k)})

    # Generic tests.
    cross_kernel_tests(k, X0, X1, X2)
    binary_map_tests(k, X0)
    pairwise_tests(k, X0; rtol=rtol)

    # Kernels should be symmetric for same arguments.
    @test pairwise(k, X0) isa LazyPDMat
    @test pairwise(k, X0) ≈ pairwise(k, X0)'

    # k(x, x′) == k(x′, x)
    @test map(k, X0, X1) ≈ map(k, X1, X0)
    @test pairwise(k, X0, X2) ≈ pairwise(k, X2, X0)'

    # Should be (approximately) positive definite.
    @test all(eigvals(Matrix(pairwise(k, X0))) .> -1e-9)

    # length should equal length of first side.
    @test length(k) == size(k, 1)
    @test size(k, 1) == size(k, 2)
end


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

"""
    ClaimDevicesByClaimCode()

Adds device(s) to your account (i.e., claim one or more devices) if and only if you
 received a claim code with the device(s).

# Required Parameters
- `claimCode`: The claim code, starting with \"C-\", as provided by the device manufacturer.

"""
claim_devices_by_claim_code(claimCode; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/claims/$(claimCode)"; aws_config=aws_config)
claim_devices_by_claim_code(claimCode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/claims/$(claimCode)", args; aws_config=aws_config)

"""
    DescribeDevice()

Given a device ID, returns a DescribeDeviceResponse object describing the
 details of the device.

# Required Parameters
- `deviceId`: The unique identifier of the device.

"""
describe_device(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)"; aws_config=aws_config)
describe_device(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)", args; aws_config=aws_config)

"""
    FinalizeDeviceClaim()

Given a device ID, finalizes the claim request for the associated device.
 Claiming a device consists of initiating a claim, then publishing a device event,
 and finalizing the claim. For a device of type button, a device event can
 be published by simply clicking the device.

# Required Parameters
- `deviceId`: The unique identifier of the device.

# Optional Parameters
- `tags`: A collection of key/value pairs defining the resource tags. For example, {
 \"tags\": {\"key1\": \"value1\", \"key2\": \"value2\"} }. For more information, see AWS
 Tagging Strategies.

"""
finalize_device_claim(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/finalize-claim"; aws_config=aws_config)
finalize_device_claim(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/finalize-claim", args; aws_config=aws_config)

"""
    GetDeviceMethods()

Given a device ID, returns the invokable methods associated with the device.

# Required Parameters
- `deviceId`: The unique identifier of the device.

"""
get_device_methods(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/methods"; aws_config=aws_config)
get_device_methods(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/methods", args; aws_config=aws_config)

"""
    InitiateDeviceClaim()

Given a device ID, initiates a claim request for the associated device.
 Claiming a device consists of initiating a claim, then publishing a device event,
 and finalizing the claim. For a device of type button, a device event can
 be published by simply clicking the device.

# Required Parameters
- `deviceId`: The unique identifier of the device.

"""
initiate_device_claim(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/initiate-claim"; aws_config=aws_config)
initiate_device_claim(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/initiate-claim", args; aws_config=aws_config)

"""
    InvokeDeviceMethod()

Given a device ID, issues a request to invoke a named device method (with possible
 parameters). See the \"Example POST\" code snippet below.

# Required Parameters
- `deviceId`: The unique identifier of the device.

# Optional Parameters
- `deviceMethod`: The device method to invoke.
- `deviceMethodParameters`: A JSON encoded string containing the device method request parameters.
"""
invoke_device_method(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/devices/$(deviceId)/methods"; aws_config=aws_config)
invoke_device_method(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/devices/$(deviceId)/methods", args; aws_config=aws_config)

"""
    ListDeviceEvents()

Using a device ID, returns a DeviceEventsResponse object containing an
 array of events for the device.

# Required Parameters
- `deviceId`: The unique identifier of the device.
- `fromTimeStamp`: The start date for the device event query, in ISO8061 format. For example,
 2018-03-28T15:45:12.880Z

- `toTimeStamp`: The end date for the device event query, in ISO8061 format. For example,
 2018-03-28T15:45:12.880Z

# Optional Parameters
- `maxResults`: The maximum number of results to return per request. If not set, a default value of
 100 is used.
- `nextToken`: The token to retrieve the next set of results.
"""
list_device_events(deviceId, fromTimeStamp, toTimeStamp; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/events", Dict{String, Any}("fromTimeStamp"=>fromTimeStamp, "toTimeStamp"=>toTimeStamp); aws_config=aws_config)
list_device_events(deviceId, fromTimeStamp, toTimeStamp, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices/$(deviceId)/events", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("fromTimeStamp"=>fromTimeStamp, "toTimeStamp"=>toTimeStamp), args)); aws_config=aws_config)

"""
    ListDevices()

Lists the 1-Click compatible devices associated with your AWS account.

# Optional Parameters
- `deviceType`: The type of the device, such as \"button\".
- `maxResults`: The maximum number of results to return per request. If not set, a default value of
 100 is used.
- `nextToken`: The token to retrieve the next set of results.
"""
list_devices(; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices"; aws_config=aws_config)
list_devices(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/devices", args; aws_config=aws_config)

"""
    ListTagsForResource()

Lists the tags associated with the specified resource ARN.

# Required Parameters
- `resource-arn`: The ARN of the resource.

"""
list_tags_for_resource(resource_arn; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/tags/$(resource-arn)"; aws_config=aws_config)
list_tags_for_resource(resource_arn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("GET", "/tags/$(resource-arn)", args; aws_config=aws_config)

"""
    TagResource()

Adds or updates the tags associated with the resource ARN. See AWS IoT 1-Click Service Limits for the maximum number of tags allowed per
 resource.

# Required Parameters
- `resource-arn`: The ARN of the resource.
- `tags`: A collection of key/value pairs defining the resource tags. For example, {
 \"tags\": {\"key1\": \"value1\", \"key2\": \"value2\"} }. For more information, see AWS
 Tagging Strategies.

"""
tag_resource(resource_arn, tags; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/tags/$(resource-arn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resource_arn, tags, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("POST", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), args)); aws_config=aws_config)

"""
    UnclaimDevice()

Disassociates a device from your AWS account using its device ID.

# Required Parameters
- `deviceId`: The unique identifier of the device.

"""
unclaim_device(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/unclaim"; aws_config=aws_config)
unclaim_device(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/unclaim", args; aws_config=aws_config)

"""
    UntagResource()

Using tag keys, deletes the tags (key/value pairs) associated with the specified
 resource ARN.

# Required Parameters
- `resource-arn`: The ARN of the resource.
- `tagKeys`: A collections of tag keys. For example, {\"key1\",\"key2\"}

"""
untag_resource(resource_arn, tagKeys; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("DELETE", "/tags/$(resource-arn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resource_arn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("DELETE", "/tags/$(resource-arn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateDeviceState()

Using a Boolean value (true or false), this operation
 enables or disables the device given a device ID.

# Required Parameters
- `deviceId`: The unique identifier of the device.

# Optional Parameters
- `enabled`: If true, the device is enabled. If false, the device is
 disabled.
"""
update_device_state(deviceId; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/state"; aws_config=aws_config)
update_device_state(deviceId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = iot_1click_devices_service("PUT", "/devices/$(deviceId)/state", args; aws_config=aws_config)


text:
#
#----------------------------------------------------------------------------#
abstract type AbstractPOU end # partition of unity
#----------------------------------------------------------------------------#
export nnPOU
#--------------------------------------#
struct nnPOU <: AbstractPOU
    nprt::Integer
    NN::Chain
end
#--------------------------------------#
function nnPOU(nprt::Integer;
               depth=2,
               act=relu,
               init=Flux.glorot_uniform)

    """
     relu activation function for sharp (piecewise linear)
     partition functions necessary for representing piecewise
     continuous functions
    """

    NN = Chain(Dense(init(nprt,1),init(nprt),act)
              ,Chain([Dense(init(nprt,nprt),init(nprt),act)
                                           for i in 1:(depth-1)]...)
              ,softmax)

    return nnPOU(nprt,NN)
end
#--------------------------------------#
Flux.@functor nnPOU
DiffEqFlux.initial_params(pou::nnPOU) = Flux.destructure(pou)[1]
#--------------------------------------#
(pou::nnPOU)(x) = pou.NN(x)
#----------------------------------------------------------------------------#
export nnPOU2D
#--------------------------------------#
struct nnPOU2D <: AbstractPOU
    nprt::Integer
    sdf
    tol::Number
    NN::Chain
end
#--------------------------------------#
function nnPOU2D(nprt::Integer,sdf
                ;tol=1f-8
                ,depth=4
                ,act=relu
                ,init=Flux.glorot_uniform)

    """
     integrate signed distance function knowledge in neural network
    """

    NN = Chain(Dense(init(nprt,2),init(nprt),act)
              ,Chain([Dense(init(nprt,nprt),init(nprt),act)
                                           for i in 1:(depth-1)]...)
              ,softmax)

    return nnPOU2D(nprt,sdf,tol,NN)
end
#--------------------------------------#
Flux.@functor nnPOU2D
DiffEqFlux.initial_params(pou::nnPOU2D) = Flux.destructure(pou)[1]
#--------------------------------------#
function (pou::nnPOU2D)(x,y)
    @unpack sdf, tol, NN = pou

    ϕ = NN(vcat(x,y))

    Ω = Zygote.ignore() do
        Ω = sdf(x,y) .<= 0f0 + tol
        return Ω
    end

    ϕ = ϕ .* Ω # clip partitions at domain boundary

    return ϕ
end
#----------------------------------------------------------------------------#
export femPOU
#--------------------------------------#
"""
 FEM partition of unity with arbitrary shape function
"""
struct femPOU <: AbstractPOU
    nprt::Integer
    shape::Function
    xe::AbstractVector # interior points
    x0::Number         # end points
    x1::Number
end

#--------------------------------------#
function femPOU(nprt::Integer
               ,shape = x -> x # linear FEM hat function
               ;x0::Number=-1f0
               ,x1::Number= 1f0)

    """
     shape: shape function defined on [0,1]
     s.t shape(0.0) = 0, shape(1.0) = 1
    """

    xg = Array(range(x0,x1,length=nprt))

    return femPOU(nprt,shape,xg[2:end-1],x0,x1)
end
#--------------------------------------#
Flux.@functor femPOU
DiffEqFlux.initial_params(pou::femPOU) = Flux.destructure(pou)[1]
#--------------------------------------#
function (pou::femPOU)(x)
    @unpack shape,nprt,xe,x0,x1 = pou

    xg = addEndPts(xe,x0,x1)

    # allow partitions to collapse onto themselves
#   xg = sort(xg)

    n  = length(xg)
    h  = diff(xg)
    hi = 1.0 ./ h

    RR,b2 = Zygote.ignore() do
        RR = restrictUpper(n) + restrictLower(n)
        RR = RR'
        b2 = zero(xg); b2[1] = 1
        return RR,b2
    end

    hh = -RR*hi
    W2 = Tdg(hi,hh,hi)

    ϕ = Chain(x -> relu.(x .- xg)
             ,u -> (W2*u .+ b2)
             ,@. shape)

    u = ϕ(x)

    return u
end
#----------------------------------------------------------------------------#
export femPOU2D
#--------------------------------------#
"""
 2D FEM partition of unity with arbitrary shape function
"""
struct femPOU2D <: AbstractPOU
    nprt::Integer
    shape::Function
    sdf  ::Function # range [-1,0]
end
#--------------------------------------#
Flux.@functor femPOU2D
DiffEqFlux.initial_params(pou::femPOU2D) = Flux.destructure(pou)[1]
#--------------------------------------#
function (pou::femPOU2D)(x,y)

    @unpack nprt, shape, sdf = pou

    N = size(x)[end]

    d  = sdf(x,y) # range [-1,0]
    dd = 1 .+ d   # range [0 ,1]
    s = shape(dd)

    u = vcat(s, 1 .- s)

    return u
end
#----------------------------------------------------------------------------#
export indicatorPOU
#--------------------------------------#
"""
 Indicator function over grid
"""
struct indicatorPOU <: AbstractPOU
    nprt::Integer
    xe::AbstractVector # interior points
    x0::Number         # end points
    x1::Number
end

function indicatorPOU(nprt
                     ;x0::Number=-1f0
                     ,x1::Number= 1f0)

    xg = Array(range(x0,x1,length=nprt+1))

    return indicatorPOU(nprt,xg[2:end-1],x0,x1)
end
#--------------------------------------#
Flux.@functor indicatorPOU
DiffEqFlux.initial_params(pou::indicatorPOU) = Flux.destructure(pou)[1]
#--------------------------------------#
function (pou::indicatorPOU)(x) # [1,N]

    @unpack nprt,xe,x0,x1 = pou

    N = size(x)[end]
    y = zeros(nprt,N)

    xg = addEndPts(xe,x0,x1) # [nprt+1]

    Ru,Rl = Zygote.ignore() do
        Ru = restrictUpper(nprt+1)
        Rl = restrictLower(nprt+1)
        return Ru,Rl
    end

    xm = Rl*xg # [nprt]
    xM = Ru*xg

    """
    julia> xp
    1×1000 adjoint(::Vector{Float32}) with eltype Float32:
     -1.0  -0.997998  -0.995996  …  0.993994  0.995996  0.997998  1.0
    julia> sign.(xp)
    1×1000 Matrix{Float32}:
     -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  …  1.0  1.0  1.0  1.0  1.0  1.0  1.0

    #########

    julia> fun(x)=sum(x .< 1);ps=Flux.params(xp);
    julia> gs = gradient(() -> fun(xp),ps);
    julia> gs[xp]
    gs[xp]
    > nothing

    #########

    julia> fun(x)=sum(sign.(x));ps=Flux.params(xp);
    julia> gs = gradient(() -> fun(xp),ps);
    julia> gs[xp]
    > 1×1000 Matrix{Float32}:
     0.0  0.0  0.0  0.0  0.0  0.0  0.0  …  0.0  0.0  0.0  0.0  0.0  0.0  0.0
    """
    y = xm .<= x .<= xM # [nprt,N]

    return y
end
#--------------------------------------#


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

using Markdown
using InteractiveUtils

# ╔═╡ 5e8c4976-3727-11eb-2801-0313fd547ac2
file = open(f->read(f, String), "day_1.txt")

# ╔═╡ bd9f511c-3727-11eb-38e3-ad2be4ac32fc
begin
	#separing each line
	input = eachmatch(r"(.+)\n",file)

	#getting first match of regex and converting to int
	input = map(x-> parse(Int64, x.captures[1]), input)
end

# ╔═╡ b3710000-372a-11eb-25f0-4fe8d0aa2162
sum = filter(x->x<2020, [i+j for i in input,j in input])

# ╔═╡ a064fb9e-3737-11eb-06e6-fb89f785aba8
needs = [(2020-x) for x in sum]

# ╔═╡ 081e453e-372d-11eb-09ab-471d8bb8881f
inter = intersect(needs,input)

# ╔═╡ cd91eaa0-3737-11eb-2b49-b9168dab6e3f
result =inter[1]*inter[2]*inter[3]

# ╔═╡ 33efb69e-3738-11eb-3c85-8df7bac1208e
1051 * 897 * 72

# ╔═╡ Cell order:
# ╠═5e8c4976-3727-11eb-2801-0313fd547ac2
# ╠═bd9f511c-3727-11eb-38e3-ad2be4ac32fc
# ╠═b3710000-372a-11eb-25f0-4fe8d0aa2162
# ╠═a064fb9e-3737-11eb-06e6-fb89f785aba8
# ╠═081e453e-372d-11eb-09ab-471d8bb8881f
# ╠═cd91eaa0-3737-11eb-2b49-b9168dab6e3f
# ╠═33efb69e-3738-11eb-3c85-8df7bac1208e


text:
module Test_pslifting

using PeriodicSystems
using DescriptorSystems
using Symbolics
using Test
using LinearAlgebra
using ApproxFun

println("Test_liftings")

# # selected sequence to illustrate failure of BlockMatrices.jl v0.16.16
# using LinearAlgebra
# using ApproxFun
# Af = Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],Fourier(0..2π));
# D = Derivative(domain(Af));
# ND = [D 0I; 0I D];
# Aop = Af - ND;
# NA = 322
# RW = Aop[1:NA,1:NA]

@testset "test_liftings" begin

# using Floquet based approach
At = PeriodicFunctionMatrix(t -> [0 1; -10*cos(t) -24-10*sin(t)],2pi);
@time ev = psceig(At; reltol = 1.e-10)
@time ev = psceig(At,500; reltol = 1.e-10)
@test sort(ev) ≈ sort([0;-24])

# using Fourier series
Afun = FourierFunctionMatrix(Fun(t -> [0 1; -10*cos(t) -24-10*sin(t)],Fourier(0..2π)))
ev1 = psceig(Afun,50)
@test isapprox(sort(real(ev)),sort(real(ev1)),rtol=1.e-6) && norm(imag(ev1)) < 1.e-10

ev3 = psceig(Afun,60)
@test isapprox(sort(real(ev)),sort(real(ev3)),rtol=1.e-6) && norm(imag(ev3)) < 1.e-10

# using Toeplitz operator truncation
Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [0 1; -10*cos(t) -24-10*sin(t)],2pi));
ev2 = psceig(Ahr,50)
@test sort(ev) ≈ sort(real(ev2[sortperm(imag(ev2),by=abs)][1:2]))

# example Zhou, Hagiwara SCL 2002 period pi/2 and pi
# using Floquet based approach
At = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi/2);
At2 = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi);
@time ev = psceig(At; reltol = 1.e-10)
@time ev = psceig(At,500; reltol = 1.e-10) # complex characteristic exponents
@test sort(real(ev)) ≈ sort([-1;-2]) && sort(imag(ev)) ≈ [2;2]
@time ev2 = psceig(At2; reltol = 1.e-10)
@time ev2 = psceig(At2,500; reltol = 1.e-10) # real characteristic exponents
@test sort(ev2) ≈ sort([-1;-2])

# using Fourier series
Afun = FourierFunctionMatrix(Fun(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],Fourier(0..π/2)))
ev3 = psceig(Afun)
@test sort(real(ev3)) ≈ sort([-1;-2]) && sort(imag(ev3)) ≈ [2;2]

# using Toeplitz operator truncation period pi/2
Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi/2));
Asym = convert(PeriodicSymbolicMatrix,Ahr); Asym.F
ev3 = psceig(Ahr,5)
@test sort(real(ev3)) ≈ sort([-1;-2]) && sort(imag(ev3)) ≈ [2;2]

# simple period
ev4 = psceig(Ahr,50)
@test sort(real(ev4)) ≈ sort([-1;-2]) && sort(imag(ev4)) ≈ [2;2]

# double period
ev5 = psceig(Ahr,50; P = 2)
@test sort(real(ev5)) ≈ sort([-1;-2]) && norm(imag(ev5)) < 1.e-10

# using Fourier series truncation period pi
ev4 = psceig(Afun,P = 2)
@test sort(real(ev4)) ≈ sort([-1;-2]) && norm(imag(ev4)) < 1.e-10

# using Toeplitz operator truncation period pi
Ahr2 = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi));
Asym2 = convert(PeriodicSymbolicMatrix,Ahr2); Asym2.F
ev3 = psceig(Ahr2,5)
@test sort(real(ev3)) ≈ sort([-1;-2]) && norm(imag(ev3)) < 1.e-10

ev4 = psceig(Ahr,50; P = 2)
@test sort(real(ev4)) ≈ sort([-1;-2]) && norm(imag(ev4)) < 1.e-10

# constant matrix case
Q = convert(HarmonicArray,PeriodicFunctionMatrix([-1 0; 0 -2],pi/2));
N = 5; ev1 = psceig(Q,N)
@test sort(real(ev1)) ≈ sort([-1;-2]) && norm(imag(ev1)) < 1.e-10

QF = convert(FourierFunctionMatrix,Q)
N = 5; ev2 = psceig(QF,N)
@test sort(real(ev2)) ≈ sort([-1;-2]) && norm(imag(ev2)) < 1.e-10

# Floquet approach
A = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi/2);
@time ev = psceig(A; reltol = 1.e-10)
@time ev = psceig(A,50; reltol = 1.e-10)

# Fourier approach
Afun = FourierFunctionMatrix(Fun(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],Fourier(0..π)))
ev1 = psceig(Afun)
@test sort(real(ev)) ≈ sort(real(ev1)) && norm(imag(ev1)) < 1.e-10

ev2 = psceig(Afun,5)
@test sort(real(ev)) ≈ sort(real(ev2)) && norm(imag(ev2)) < 1.e-10

ev3 = psceig(Afun,20)
@test sort(real(ev)) ≈ sort(real(ev3)) && norm(imag(ev3)) < 1.e-10

# # lossy Mathieu differential equation
# k = -1.; ξ = 0.05; β =0.2; ωh = 2; # unstable feedback
# k = -.5; ξ = 0.05; β =0.2; ωh = 2; # stable feedback
# Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [0 1; k*(1-2β*cos(ωh*t)) -2*ξ],2pi));
# Asym = convert(PeriodicSymbolicMatrix,Ahr); Asym.F
# N = 5; ev = eigvals(hr2bt(Ahr,N)-phasemat(Ahr,N))
# ev[sortperm(imag(ev),by=abs)][1:2]

# Bitanti-Colaneri's book p.26
a = PeriodicFunctionMatrix(t -> [-1+sin(t) 0; 1-cos(t) -3],2*pi);
b = PeriodicFunctionMatrix(t ->  [-1-cos(t);2-sin(t)], 2*pi);
c = [0 1]; d = [0];

# using Harmonic Array based lifting
psyschr = ps(HarmonicArray,a,b,c,d)
ev = psceig(psyschr.A,20)

syshr = ps2fls(psyschr,20)
p = gpole(syshr,atol=1.e-7);  p = p[sortperm(imag(p),by=abs)][1:2];
@test sort(real(p)) ≈ sort(real(ev))  && norm(imag(p)) < 1.e-10

z = gzero(syshr,atol=1.e-7); z = z[isfinite.(z)]; # Question: How to handle infinite zeros?
ζ = z[sortperm(imag(z),by=abs)][1]

Z = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [(-2+3*sin(t))/(2-sin(t))],2*pi))
ρ = real(Z.values[:,:,1])[1,1]
@test real(ζ) ≈ ρ

# using Fourier Function Matrix based lifting
psysc = ps(FourierFunctionMatrix,a,b,c,d);
ev1 = psceig(psysc.A,30)

sys = ps2frls(psysc,20);
p = gpole(sys,atol=1.e-7);  p = p[sortperm(imag(p),by=abs)][1:4];
@test minimum(abs.(p.-ev1[1])) < 1.e-10 && minimum(abs.(p.-ev1[2])) < 1.e-10

z = gzero(sys,atol=1.e-7); z = z[isfinite.(z)]; # Question: How to handle infinite zeros?
z = z[sortperm(imag(z),by=abs)][1:4];
@test minimum(abs.(z.-ρ)) < 1.e-10

# Zhou-Hagiwara Automatica 2002
β = 0.5
a1 = PeriodicFunctionMatrix(t -> [-1-sin(2*t)^2 2-0.5*sin(4*t); -2-0.5*sin(4*t) -1-cos(2*t)^2],pi);
#γ(t) = mod(t,pi) < pi/2 ? sin(2*t) : 0
γ = chop(Fun(t -> mod(t,pi) < pi/2 ? sin(2*t) : 0, Fourier(0..pi)),1.e-7);
b1 = PeriodicFunctionMatrix(t ->  [0; 1-2*β*(mod(t,float(pi)) < pi/2 ? sin(2*t) : 0 )], pi);
b1 = PeriodicFunctionMatrix(t ->  [0; 1-2*β*γ(t)], pi);
c = [1 1]; d = [0];

# using Harmonic Array based lifting
psyschr = ps(HarmonicArray,a1,b1,c,d);
ev = psceig(psyschr.A)

syshr = ps2fls(psyschr,60);
p = gpole(syshr,atol=1.e-7);  p = p[sortperm(imag(p),by=abs)][1:2]
@test sort(real(p)) ≈ sort(real(ev))  && norm(imag(p)) < 1.e-10

z = gzero(syshr,atol=1.e-7); z = z[isfinite.(z)]; #  Question: How to handle infinite zeros?
ζ = z[sortperm(imag(z),by=abs)][1]
@test ζ ≈ -3.5

# using Fourier Function Matrix based lifting
psysc = ps(FourierFunctionMatrix,a1,b1,c,d);
ev = psceig(psysc.A)

sys = ps2frls(psysc,60);
p = gpole(sys,atol=1.e-7);  p = p[sortperm(imag(p),by=abs)][1:2]
@test sort(real(p)) ≈ sort(real(ev))  && norm(imag(p)) < 1.e-10

z = gzero(sys,atol=1.e-7); z = z[isfinite.(z)]; #  Question: How to handle infinite zeros?
z = z[sortperm(imag(z),by=abs)][1]
@test z[1] ≈ -3.5

# # Ziegler's column

# β = 2pi; η =0.5; λ = 0.5;
# M = [1 3/8; 3/2 1];
# K0 = [3/8 3/16; -3/4 3/4]; K1 = λ*[-1 η; 0 4η-4];
# A = PeriodicFunctionMatrix(τ -> [eye(2) zeros(2,2); zeros(2,2) -M\(K0+cos(β*τ).*K1)],2*pi/β)

# # Mathieu equation: determination of minimum eigenvalue
# q = 2; a = 2; a = -1.513956885056448; a = -1
# Ahr = convert(HarmonicArray,PeriodicFunctionMatrix(t -> [0 1; a-2*q*cos(2*t) 0],pi));
# Asym = convert(PeriodicSymbolicMatrix,Ahr); Asym.F
# ev = psceig(Ahr)
# N = 20; ev = eigvals(hr2bt(Ahr,N)-phasemat(Ahr,N))
# ev[sortperm(imag(ev),by=abs)][1:2]

# N = 5; ev1 = eigvals(hr2bt1(Ahr,N,2)-phasemat1(Ahr,N,2))
# ev1[sortperm(imag(ev1),by=abs)][1:2]

# discrete liftings
Ad = PeriodicMatrix([[1. 0; 0 0], [1 1;1 1], [0 1; 1 0]], 6, nperiod = 2);
Bd = PeriodicMatrix( [[ 1; 0 ], [ 1; 1]] ,2);
Cd = PeriodicMatrix( [[ 1 1 ], [ 1 0]] ,2);
Dd = PeriodicMatrix( [[ 1 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
sys = ps2ls(psys)
sys1 = ps2ls(psys, 7, ss = true)
@test iszero(sys-sys1)

Ad = PeriodicMatrix([[1. 0], [1;1]],2);
Bd = PeriodicMatrix( [[ 1 ], [ 1; 1]] ,2);
Cd = PeriodicMatrix( [[ 1 1 ], [ 1 ]] ,2);
Dd = PeriodicMatrix( [[ 1 ]], 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
sys = ps2ls(psys)
sys1 = ps2ls(psys, 7, ss = true)
@test iszero(sys-sys1)

sys = ps2ls(psys, 2)
sys1 = ps2ls(psys, 8, ss = true)
@test iszero(sys-sys1)

sys = ps2ls(psys, cyclic = true)
sys1 = ps2ls(psys, 3, cyclic = true)
@test iszero(sys-sys1)

Ad = PeriodicArray(rand(Float32,2,2,10),10);
Bd = PeriodicArray(rand(2,1,2),2);
Cd = PeriodicArray(rand(1,2,3),3);
Dd = PeriodicArray(rand(1,1,1), 1);
psys = PeriodicStateSpace(Ad,Bd,Cd,Dd);
sys = ps2ls(psys)
sys1 = ps2ls(psys, ss = true)
@test iszero(sys-sys1)

sys = ps2ls(psys, 2)
sys1 = ps2ls(psys, 2, ss = true)
@test iszero(sys-sys1)

sys = ps2ls(psys, cyclic = true)
sys1 = ps2ls(psys, 31, cyclic = true)
@test iszero(sys-sys1)

end # test

end # module


text:
import Base: exp, log
import Base: .+, +, .-, -, .*, *
import Base: broadcast, broadcast!

macro elemwise(op)
    quote
        op = $op
        y = CuArray{T}(size(x))
        t = ctype(T)
        f = @nvrtc """
        $array_h
        __global__ void f(Array<$t,$N> x, Array<$t,$N> y) {
            int idx = blockIdx.x * blockDim.x + threadIdx.x;
            if (idx < y.length()) {
                y(idx) = $op(x(idx));
            }
        } """
        f(x, y, dx=length(y))
        y
    end
end

function exp{T,N}(x::CuArray{T,N})
    f = @elemwise (x,y) """
    y(idx) = exp(x(idx));
    """
    f(x)
end

function +{T,N}(x1::CuArray{T,N}, x2::CuArray{T,N})
    y = similar(x1)
    f = @elemwise (x1,x2,y) """
    y(idx) = x1(idx) + x2(idx);
    """
    f(x)
end

for op in (:exp, :log)
    @eval begin
        function $op{T,N}(x::AbstractCuArray{T,N})
            op = $op
            y = CuArray{T}(size(x))
            t = ctype(T)
            f = @nvrtc """
            $array_h
            __global__ void f(Array<$t,$N> x, Array<$t,$N> y) {
                int idx = blockIdx.x * blockDim.x + threadIdx.x;
                if (idx < y.length()) {
                    y(idx) = $op(x(idx));
                }
            } """
            f(x, y, dx=length(y))
            y
        end
    end
end

for op in (:+, :-)
    @eval begin
        function $op{T,N}(x1::AbstractCuArray{T,N}, x2::AbstractCuArray{T,N})
            size(x1) == size(x2) || throw(DimensionMismatch())
            op = $op
            y = CuArray{T}(size(x1))
            t = ctype(T)
            f = @nvrtc """
            $array_h
            __global__ void f(Array<$t,$N> x1, Array<$t,$N> x2, Array<$t,$N> y) {
                int idx = blockIdx.x * blockDim.x + threadIdx.x;
                if (idx < y.length()) {
                    y(idx) = x1(idx) $op x2(idx);
                }
            } """
            f(x1, x2, y, dx=length(y))
            y
        end
    end
end

for op in (:+, :-, :*)
    @eval begin
        function broadcast!{T,N}(::typeof($op), y::AbstractCuArray{T,N},
            x1::AbstractCuArray{T,N}, x2::AbstractCuArray{T,N})
            op = $op
            t = ctype(T)
            f = @nvrtc """
            $array_h
            __global__ void f(Array<$t,$N> y, Array<$t,$N> x1, Array<$t,$N> x2) {
                int idx = blockIdx.x * blockDim.x + threadIdx.x;
                if (idx < y.length()) {
                    int subs[$N];
                    y.idx2sub(idx, subs);
                    y(subs) = x1(subs) $op x2(subs);
                }
            } """
            f(y, x1, x2, dx=length(y))
            y
        end
    end
end

for (op1,op2) in ((:.+,:+), (:.-,:-))
    @eval begin
        function $op1{T,N}(x1::AbstractCuArray{T,N}, x2::AbstractCuArray{T,N})
            dims = ntuple(i -> max(size(x1,i),size(x2,i)), N)
            y = CuArray{T}(dims)
            broadcast!($op2, y, x1, x2)
        end
    end
end

*(x1::CuMatrix, x2::CuMatrix) = BLAS.gemm(x1, x2)


text:
mutable struct Post
    opinion::Float64
    weight::Float64
    source_agent::Int64
    published_at::Int64
    seen_by::Array{Int64, 1}
    function Post(opinion, weight, source_agent, published_at)
        # check if opinion value is valid
        if opinion < -1 || opinion > 1
            error("invalid opinion value")
        end
        if weight < 0
            error("invalid weight value")
        end
        new(
            opinion,
            weight,
            source_agent,
            published_at,
            Int64[]
        )
    end
end

Base.:<(x::Post, y::Post) = x.weight < y.weight

# suppress output of include()
;


text:
using indentationTools
using Test

include("testFunctions.jl")

#importNI_forceDisplacementData("test/testData/anonNI.TXT")
@testset "importNI_forceDisplacementData " begin
    @test begin
        2000 == length(importNI_forceDisplacementData("test/testData/anonNI.TXT"))
    end
end

@testset "Test readIBW                   " begin
    # TO DO:
    # ERROR HANDLING
    @test  begin # "Read back written data: small array"
        filename = "Test1.ibw"
        inputData = Float32.(rand(100,3))
        writeIBW(filename,inputData)
        all(reshape(inputData, : , 1) .== readIBW(filename))
    end
    rm("Test1.ibw")
    @test  begin # "Read back written data: big array"
        filename = "Test2.ibw"
        inputData = Float32.(rand(10_000,3))
        writeIBW(filename,inputData)
        all(reshape(inputData, : , 1) .== readIBW(filename))
    end
    rm("Test2.ibw")
    @test begin # "Data supplied as a Float64 matrix"
        filename = "Test3.ibw"
        inputData = Float64.(rand(100,3))
        writeIBW(filename,inputData)
        all(Float32.(reshape(inputData, : , 1)) .== readIBW(filename))
    end
    rm("Test3.ibw")

    ### Here comes a series of regression tests designed on real data files.
    ### Obviously any failing test means only that behavior has changed, not
    ### necessarily that the previous implementation was correct. This is
    ### because the real .IBW file format is "black-box" to me.
    @test begin
        # Regression test 1
        inputVec = readIBW("test/testData/anonMG.ibw")
        length(inputVec) == 254541
    end
    @test begin
        # Regression test 2
        inputVec = readIBW("test/testData/anonMG.ibw")
        inputVec[12345] == -3.404614f-6
    end
    @test begin
        # Regression test 3
        inputVec = readIBW("test/testData/anonZF.ibw")
        length(inputVec) == 254466
    end
    @test begin
        # Regression test 4
        inputVec = readIBW("test/testData/anonZF.ibw")
        inputVec[end] == -2.6649975f-6
    end
    @test begin
        # Regression test 5
        inputVec = readIBW("test/testData/anonPMMA.ibw")
        length(inputVec) == 254961
    end
    @test begin
        # Regression test 6
        inputVec = readIBW("test/testData/anonPMMA.ibw")
        inputVec[1] == -2.3641041f-6
    end
    @test begin
        # Regression test 7
        inputVec = readIBW("test/testData/anonCell.ibw")
        length(inputVec) == 4299111
    end
    @test begin
        # Regression test 8
        inputVec = readIBW("test/testData/anonCell.ibw")
        inputVec[1337] == 1.5439662f-7
    end
end

@testset "Test IBWtoTXT                  " begin
    @test begin
        filename = "Test1.ibw"
        inputData = Float32.([collect(1:100) collect(101:200) collect(201:300)])
        writeIBW(filename,inputData)
        all(inputData[:,3:-1:2] .== IBWtoTXT(filename))
    end
    rm("Test1.ibw")
    @test begin
        filename = "Test1.ibw"
        inputData = Float32.(rand(100,3))
        writeIBW(filename,inputData)
        all(inputData[:,3:-1:2] .== IBWtoTXT(filename))
    end
    rm("Test1.ibw")
end

@testset "Test offsetAndDriftCompensation" begin
    ############################################
    # First data set:
    #   - No noise
    #   - Linear interpolation between keypoints
    #   -
    sampleRate = 1000.0 # Hz
    #                 TIME   FORCE    DISPLACEMENT
        matOfKeyPoints = [ 0.0    100.0    -200.0        # Starting position
                           5.0    100.0     100.0        # Position of contact
                           6.0   1500.0     110.0        # End of ramp / Start of hold
                          16.0   1500.0     125.0        # End of hold / Start of unload
                          17.0    175.0     105.0        # End of unload / Start of thermal hold
                          47.0    175.0     101.0        # End of thermal hold / Start of retraction
                          48.0    100.0     100.0        # End of retraction / Start of exit
                          50.0    100.0    -100.0]       # End of exit & End of experiment
    xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate)
    xy , basefit , zeroLine , rampStartIdx, xZero = offsetAndDriftCompensation(Float32.(xy))
    @test begin
        isapprox(zeroLine[1] , (matOfKeyPoints[2,2]-matOfKeyPoints[1,2])/(matOfKeyPoints[2,3]-matOfKeyPoints[1,3]), atol = 1e-6)
    end
    @test begin
        isapprox(zeroLine[2] , 0.5*(matOfKeyPoints[1,2]+matOfKeyPoints[2,2]), atol = 1e-6)
    end
    @test begin
        isapprox(rampStartIdx , sampleRate.*matOfKeyPoints[2,1], rtol = 0.01)
    end
    @test begin
        isapprox(Float32(xZero), matOfKeyPoints[2,3], rtol = 0.01)
    end
    ############################################
    # Second data set:
    #   - Noise
    #   - Linear interpolation between keypoints
    #   -
    sampleRate = 1000.0 # Hz
    #                 TIME   FORCE    DISPLACEMENT
        matOfKeyPoints = [ 0.0    100.0    -200.0        # Starting position
                           5.0    100.0     100.0        # Position of contact
                           6.0   1500.0     110.0        # End of ramp / Start of hold
                          16.0   1500.0     125.0        # End of hold / Start of unload
                          17.0    175.0     105.0        # End of unload / Start of thermal hold
                          47.0    175.0     101.0        # End of thermal hold / Start of retraction
                          48.0    100.0     100.0        # End of retraction / Start of exit
                          50.0    100.0    -100.0]       # End of exit & End of experiment
    # xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate, stdOfNoise = 2.0,sizeOfNoise =[0.0 0.4 0.4])
    # xy , basefit , zeroLine , rampStartIdx, xZero = offsetAndDriftCompensation(Float32.(xy))
    # @test begin
    #     isapprox(zeroLine[1] , (matOfKeyPoints[2,2]-matOfKeyPoints[1,2])/(matOfKeyPoints[2,3]-matOfKeyPoints[1,3]), atol = 0.02)
    # end
    # @test begin
    #     isapprox(zeroLine[2] , 0.5*(matOfKeyPoints[1,2]+matOfKeyPoints[2,2]), atol = 10.0)
    # end
    # @test begin
    #     #println(rampStartIdx)
    #     #println(sampleRate.*matOfKeyPoints[2,1])
    #     isapprox(rampStartIdx, sampleRate.*matOfKeyPoints[2,1], atol = 25)
    # end
    # @test begin
    #     isapprox(Float32(xZero), matOfKeyPoints[2,3], rtol = 0.05)
    # end

end

@testset "Test findStartOfHold           " begin
    sampleRate = 1000.0 # Hz
    #                 TIME   FORCE    DISPLACEMENT
    matOfKeyPoints  = [ 0.0    100.0    -200.0        # Starting position
                        5.0    100.0     100.0        # Position of contact
                        6.0   1500.0     110.0        # End of ramp / Start of hold
                        16.0   1500.0     125.0        # End of hold / Start of unload
                        17.0    175.0     105.0        # End of unload / Start of thermal hold
                        47.0    175.0     101.0        # End of thermal hold / Start of retraction
                        48.0    100.0     100.0        # End of retraction / Start of exit
                        50.0    100.0    -100.0]       # End of exit & End of experiment
    xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate)
    xy , ~ , ~ , ~, ~ = offsetAndDriftCompensation(Float32.(xy))
    @test begin
        isapprox(findStartOfHold(xy, "first") , sampleRate*matOfKeyPoints[3,1], atol = 10)
    end
    @test begin
        isapprox(findStartOfHold(xy, "last") , sampleRate*matOfKeyPoints[4,1], atol = 10)
    end
    xy = generateSyntheticForceDisplacementSignal(matOfKeyPoints,sampleRate, stdOfNoise = 2.0, sizeOfNoise =[0.0 0.4 0.4])
    xy , ~ , ~ , ~, ~ = offsetAndDriftCompensation(Float32.(xy))
    @test begin
        isapprox(findStartOfHold(xy, "first") , sampleRate*matOfKeyPoints[3,1], atol = 30)
    end
    @test begin
        isapprox(findStartOfHold(xy, "last") , sampleRate*matOfKeyPoints[4,1], atol = 30)
    end
end

@testset "determineThermalCreep          " begin
    thermalHoldTime = 10000
    sampleRate = 1000
    holdTimeVals = collect(1:10000) ./ sampleRate
    thermalCreepGenerator(x,y,z) = x .+ y.*holdTimeVals.^z
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(44.0,4.0,1.1)

        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , creepVal , rtol = 0.10)
    end
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(123.0,4.0,1.1)
        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , creepVal,rtol = 0.10)
    end
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(123.0,0.20,1.1)
        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-0.2 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , creepVal,rtol = 0.10)
    end
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(123.0,4.00,0.8)
        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , creepVal,rtol = 0.15)
    end
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(0.0,4.00,0.8)
        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , creepVal,rtol = 0.15)
    end
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(123.0,0.00,0.8)
        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-0.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , creepVal,atol = 0.05)
    end
    @test begin
        xyStart = [collect( 1000.0:-1:500 ) ; 500.0.+10.0*rand(8998) ; collect(500:-1:0.0)]
        xy_hold = thermalCreepGenerator(123.0,4.00,0.0)
        creepVal , thermalHoldStartIdx = determineThermalCreep(Float32.([xy_hold[end:-1:1] xyStart]),sampleRate, thermalHoldTime, control(false, false, false), 5.0)
        isapprox(-4.0 * 0.0*holdTimeVals[end]^(0.0 - 1.0) , creepVal, atol = 0.25)
    end
end

@testset "determineCreepDuringHold       " begin
    sampleRate = 1000
    holdTimeVals = collect(1:10000) ./ sampleRate
    holdCreepGenerator(x,y,z) = x .+ y.*holdTimeVals.^z
    @test begin
        xy_hold = holdCreepGenerator(44.0,4.0,1.1)
        isapprox(4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,4.0,1.1)
        isapprox(4.0 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,0.20,1.1)
        isapprox(0.2 * 1.1*holdTimeVals[end]^(1.1 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,4.00,0.8)
        isapprox(4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(0.0,4.00,0.8)
        isapprox(4.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),rtol = 0.02)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,0.00,0.8)
        isapprox(0.0 * 0.8*holdTimeVals[end]^(0.8 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),atol = 0.01)
    end
    @test begin
        xy_hold = holdCreepGenerator(123.0,4.00,0.0)
        isapprox(4.0 * 0.0*holdTimeVals[end]^(0.0 - 1.0) , determineCreepDuringHold(xy_hold,sampleRate),atol = 0.01)
    end
end

@testset "Test hyperParameters struct    " begin
    # SamplingRate must be a positive integer.
    @test_throws DomainError hyperParameters( 0, 1200, "AP-OP", false, false, false, 0.0, 5.0 , "force")
    # SamplingRate must be a positive integer
    @test_throws DomainError hyperParameters( -2000, 1200, "AP-OP", false, false, false, 0.0, 5.0 , "force")

    # unloadingFitRange > 1.
    @test_throws DomainError hyperParameters( 1, 0, "AP-OP", false, false, false, 0.0, 5.0 , "force")
    @test_throws DomainError hyperParameters( 1, -1, "AP-OP", false, false, false, 0.0, 5.0 , "force")

    # unloadingFitFunction must be implemented.
    @test_throws DomainError hyperParameters( 1, 1200, "XP", false, false, false, 0.0, 5.0 , "force")

    # machineCompliance must not be negative.
    @test_throws DomainError hyperParameters( 1, 1200, "AP-OP", false, false, false, -0.1, 5.0 , "force")

    # noiseMultiplier must be greater than 0.
    @test_throws DomainError hyperParameters( 1, 1200, "AP-OP", false, false, false, 0.1, -5.0 , "force")
    @test_throws DomainError hyperParameters( 1, 1200, "AP-OP", false, false, false, 0.1, 0.0 , "force")

    # controlLoop must be implemented.
    @test_throws DomainError hyperParameters( 1, 1200, "AP-OP", false, false, false, 0.1, 5.0 , "forceX")
end

@testset "modulusfitter                  " begin
    ctrl = control(true, false, false)
    indentSet = metaInfoExperimentalSeries("Test2",
                                            50.0,
                                        "hemisphere",
                                        "X",
                                        217.51,
                                        "test/testData/sampleAreaFunction.txt",
                                        "test/testData/",
                                        50000,
                                        "afm")

    # Regression tests
    @test begin
        hp =  hyperParameters( 2000, 1400, "Oliver-Pharr", false, 0, 0, 0.0 , 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonMG.ibw")
        #println("Er is $Er GPa")
        isapprox(1.8696, Er ,rtol = 0.05)
    end
    @test begin
        hp =  hyperParameters( 2000, 1400, "Oliver-Pharr", true, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonMG.ibw")
        #println("Er is $Er GPa")
        isapprox(1.6, Er,rtol = 0.05)
    end
    @test begin
        hp =  hyperParameters( 2000, 1400, "Feng", false, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonMG.ibw")
        #println("Er is $Er GPa")
        isapprox( 1.9289, Er,rtol = 0.05)
    end
    @test begin
        hp = hyperParameters( 2000, 1400, "Feng", true, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonMG.ibw")
        #println("Er is $Er GPa")
        isapprox( 1.622, Er, rtol = 0.05)
    end

    indentSet = metaInfoExperimentalSeries("Test2",
                                            50.0,
                                        "hemisphere",
                                        "X",
                                        217.51,
                                        "vickers",
                                        "test/testData/",
                                        0,
                                        "ni")
     @test begin
        hp =  hyperParameters( 10, 200, "Oliver-Pharr", false, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonNI.TXT")
        #println("Er is $Er GPa")
        isapprox( 2.871314842239161, Er, rtol = 1.15)
     end
     @test begin
        hp =  hyperParameters( 10, 200, "AP-OP", false, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonNI.TXT")
        #println("Er is $Er GPa")
        isapprox( 2.9104987116718184, Er, rtol = 1.15)
     end
     @test begin
        hp =  hyperParameters( 10, 200, "Feng", false, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonNI.TXT")
        #println("Er is $Er GPa")
        isapprox( 2.8041684985519852, Er, rtol = 1.15)
    end
    @test begin
        hp =  hyperParameters( 10, 200, "Feng", true, 0, 0, 0.0, 2.0 , "force")
        Er , ~ , ~ , ~ , ~ , ~  = modulusfitter( indentSet, hp, ctrl, "anonNI.TXT")
        #println("Er is $Er GPa")
        isapprox( 2.5278698699501563 , Er, rtol = 1.15)

    end
end


text: # refer to benchmarks/run-benchmarks.jl


text:
function BlockMeans(X::Array{Float64}; alpha::Float64 = 1/2, block_size::Int = -1)
  # estimate the asymptotic variance for a stationary time series using block means
  # see v.d. Vaart, Time Series, Section 5.1
  # alpha determines block size through l = [n^alpha]
  # alpha is overruled by block_size
  if ndims(X) == 1
    n = length(X)
    X = reshape(X,(1,n))
  end
  (dimension,n) = size(X)
  if block_size < 0
    block_size = Int(floor(n^alpha))
  end
  Z = Array{Float64,2}(undef,dimension,n-block_size+1) # block means
  for i = 1:(n-block_size+1)
    Z[:,i] = mean(X[:,i:(i+block_size-1)], dims = 2)
  end
  mu = mean(X,dims=2)
  v = vec(block_size/(n-block_size+1) * sum((Z-repeat(mu,1, n-block_size+1)).^2, dims=2))
  #F(x) = 1/(n-l+1) * sum(sqrt(l)*(Z-mu).<=x)
  # return (v,F)
end

function BatchMeans(X::Matrix{Float64}; alpha::Float64 = 1/2, n_batches::Int = -1, true_mean::Vector{Float64} = Vector{Float64}(undef,0))
  # estimate the asymptotic variance for a stationary time series using batch means

  # if ndims(X) == 1
  #   n = length(X)
  #   X = reshape(X,(1,n))
  # end
  (dimension,n) = size(X)
  if (n_batches <= 0)
    batch_size = Int(floor(n^alpha))
    n_batches = Int(floor(n / n^alpha))
  else
    batch_size = Int(floor(n/n_batches))
  end
  if n_batches > n/2
    return Inf
  end
  Z = Array{Float64,2}(undef,dimension,n_batches) # batch means
  unused = n - batch_size * n_batches
  for i=1:n_batches
    Z[:,i] = mean(X[:,(unused + (i-1)*batch_size + 1):(unused + i * batch_size)], dims=2)
  end
  if length(true_mean) == 0
    vec(var(Z, dims=2)*batch_size)
  else
    sum((Z - repeat(true_mean, 1, n_batches)).^2, dims = 2)/n_batches * batch_size
  end
end

function ESS(X::AbstractArray{Float64};alpha::Float64=1/2, true_variance::Vector{Float64} = Vector{Float64}(undef,0), true_mean::Vector{Float64} = Vector{Float64}(undef,0), n_batches::Int = -1)

  # if ndims(X) == 1
  #   n = length(X)
  #   X = reshape(X,(1,n))
  # end
  asvar = vec(BatchMeans(X,alpha=alpha,n_batches=n_batches,true_mean=true_mean))
  if length(true_variance) == 0
    variance = vec(var(X,dims=2))
  else
    variance = true_variance
  end
  return size(X)[2] * variance ./ asvar
end

function SquaredJumpingDistance(X::AbstractArray{Float64})

  if ndims(X) == 1
    n_samples = length(X)
    differences = X[2:end] - X[1:end-1]
  else
    n_samples = size(X)[2]
    differences = X[:,2:end] - X[:,1:end-1]
  end
  return sumabs2(differences)/(n_samples-1)

end

function PooledVariance(AsVar::Vector{Float64}, n_samples::Vector{Int}, SampleVariance::Vector{Float64})

  return sum(n_samples .* SampleVariance ./ AsVar)/sum(n_samples ./AsVar)

end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge