text:
"""
    testpattern() -> Matrix{RGBA{N0f8}}

Load and return the provided 300x400 test image.

The returned image was specifically designed to be informative
about the effects of the applied augmentation operations. It is
thus well suited to prototype an augmentation pipeline, because it
makes it easy to see what kind of effects one can achieve with it.
"""
testpattern() = load(joinpath(@__DIR__, "..", "resources", "testpattern.png"))

function use_testpattern()
    @info("No custom image specifed. Using \"testpattern()\" for demonstration.")
    testpattern()
end

# --------------------------------------------------------------------
"""
    contiguous(A::AbstractArray)
    contiguous(A::Tuple)

Return a memory contiguous array for better performance.

Data copy only happens when necessary. For example, views returned by `view`,
`permuteddimsview` are such cases.

See also: [`plain_array`](@ref), [`plain_axes`](@ref)
"""
@inline contiguous(A::OffsetArray) = A
@inline contiguous(A::Array) = A
@inline contiguous(A::SArray) = A
@inline contiguous(A::MArray) = A
@inline contiguous(A::AbstractArray) = match_idx(collect(A), axes(A))
@inline contiguous(A::Tuple) = map(contiguous, A)

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

@inline _plain_array(A::OffsetArray) = _plain_array(parent(A))
@inline _plain_array(A::Array) = A
@inline _plain_array(A::SArray) = A
@inline _plain_array(A::MArray) = A
@inline _plain_array(A::AbstractArray) = collect(A)
@inline _plain_array(A::Tuple) = map(_plain_array, A)

"""
    plain_array(A::AbstractArray)
    plain_array(A::Tuple)

Return a memory contiguous plain array for better performance.
A plain array is either an `Array` or a `StaticArray`.

Data copy only happens when necessary. For example, views returned by `view`,
`permuteddimsview` are such cases.

See also: [`contiguous`](@ref), [`plain_axes`](@ref)
"""
@inline plain_array(A) = _plain_array(contiguous(A)) # avoid recursion

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

"""
    plain_axes(A::AbstractArray)

Generate a 1-based array from `A` without data copy.

See also: [`contiguous`](@ref), [`plain_array`](@ref)
"""
@inline plain_axes(A::Array) = A
@inline plain_axes(A::OffsetArray) = parent(A)
@inline plain_axes(A::AbstractArray) = _plain_axes(A, axes(A))
@inline plain_axes(A::SubArray) = _plain_axes(A, A.indices)

@inline function _plain_axes(A::AbstractArray{T,N}, ids::NTuple{N,Base.OneTo}) where {T, N}
    A
end

@inline function _plain_axes(A::AbstractArray, ids::Tuple{Vararg{Any}})
    view(A, axes(A)...)
end

@inline function _plain_axes(A::AbstractArray{T,N}, ids::NTuple{N, IdentityUnitRange}) where {T, N}
    view(A, map(i->i.indices, ids)...)
end

@inline function _plain_axes(A::SubArray{T,N}, ids::NTuple{N,IdentityRange}) where {T, N}
    view(parent(A), axes(A)...)
end

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

@inline match_idx(buffer::AbstractArray, inds::Tuple) = buffer
@inline match_idx(buffer::Union{Array,SubArray}, inds::NTuple{N,OffsetRange}) where {N} =
    OffsetArray(buffer, inds)

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

function indirect_axes(::Tuple{}, ::Tuple{})
    throw(MethodError(indirect_axes, ((),())))
end

@inline function indirect_axes(O::NTuple{N,Base.OneTo}, I::NTuple{N,AbstractUnitRange}) where N
    map(IdentityRange, I)
end

@inline function indirect_axes(O::NTuple{N,Base.OneTo}, I::NTuple{N,StepRange}) where N
    I
end

function indirect_axes(O::NTuple{N,AbstractUnitRange}, I::NTuple{N,AbstractUnitRange}) where N
    map((i1,i2) -> IdentityRange(UnitRange(i1)[i2]), O, I)
end

function indirect_axes(O::NTuple{N,AbstractUnitRange}, I::NTuple{N,StepRange}) where N
    map((i1,i2) -> UnitRange(i1)[i2], O, I)
end

function indirect_axes(O::NTuple{N,StepRange}, I::NTuple{N,AbstractRange}) where N
    map((i1,i2) -> i1[i2], O, I)
end

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

function indirect_view(A::AbstractArray, I::Tuple)
    view(A, indirect_axes(axes(A), I)...)
end

function indirect_view(A::SubArray{T,N,TA,<:NTuple{N,AbstractRange}}, I::Tuple) where {T,N,TA}
    view(parent(A), indirect_axes(A.indices, I)...)
end

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

function direct_axes(::Tuple{}, ::Tuple{})
    throw(MethodError(direct_axes, ((),())))
end

# TODO: Figure out why this method exists
function direct_axes(O::NTuple{N,IdentityRange}, I::NTuple{N,StepRange}) where N
    throw(MethodError(direct_axes, (O, I)))
end

@inline function direct_axes(O::NTuple{N,AbstractRange}, I::NTuple{N,AbstractUnitRange}) where N
    map(IdentityRange, I)
end

@inline function direct_axes(O::NTuple{N,AbstractRange}, I::NTuple{N,StepRange}) where N
    I
end

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

function direct_view(A::AbstractArray{T,N}, I::NTuple{N,AbstractRange}) where {T,N}
    view(A, direct_axes(axes(A), I)...)
end

function direct_view(A::SubArray{T,N,TA,<:NTuple{N,AbstractRange}}, I::NTuple{N,AbstractRange}) where {T,N,TA}
    view(A, direct_axes(A.indices, I)...)
end

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

@inline vectorize(A::AbstractVector) = A
@inline vectorize(A::Real) = A:A

@inline round_if_float(num::Integer, d) = num
round_if_float(num::AbstractFloat, d) = round(num, digits=d)
round_if_float(nums::Tuple, d) = map(num->round_if_float(num,d), nums)

function unionrange(i1::AbstractUnitRange, i2::AbstractUnitRange)
    map(min, first(i1), first(i2)):map(max, last(i1), last(i2))
end

@inline _showcolor(io::IO, T::Type{<:Number}) = print(io, T)
@inline _showcolor(io::IO, T) = ColorTypes.colorant_string_with_eltype(io, T)

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

function _2dborder!(A::AbstractArray{T,3}, val::T) where T
    ndims, h, w = size(A)
    @inbounds for j = (1,w), i = 1:h
        for d = 1:ndims
            A[d,i,j] = val
        end
    end
    @inbounds for j = 1:w, i = (1,h)
        for d = 1:ndims
            A[d,i,j] = val
        end
    end
    A
end


text:
"""
Renko chart patterns

# Methods

```
function renko(x::Array{T}; box_size::T=10.0)::Array{Int} where {T<:Real}

function renko(hlc::Matrix{T}; box_size::T=10.0, use_atr::Bool=false, n::Int=14)::Array{Int} where {T<:Real}
```

- Traditional (Constant Box Size): `renko(x::Array{T}; box_size::T=10.0)::Array{Int}`
- ATR Dynamic Box Size: `renko(hlc::Matrix{T}; box_size::T=10.0, use_atr::Bool=false, n::Int=14)::Array{Int}`

# Output

`Array{Int}` object of size Nx1 (where N is the number rows in `x`) where each element gives the Renko bar number of the corresponding row in `x`.

"""
function renko(x::Array{T}; box_size::T=10.0)::Array{Int} where {T<:Real}
    # Renko chart bar identification with traditional methodology (constant box size)
    @assert box_size != 0 "Argument `box_size` must be nonzero."
    if box_size < 0.0
        box_size = abs(box_size)
    end
    bar_id = ones(Int, size(x,1))
    ref_pt = x[1]
    @inbounds for i in 2:size(x,1)
        if abs(x[i]-ref_pt) >= box_size
            ref_pt = x[i]
            bar_id[i:end] .+= 1
        end
    end
    return bar_id
end

function renko(hlc::Matrix{T}; box_size::T=10.0, use_atr::Bool=false, n::Int=14)::Array{Int} where {T<:Real}
    # Renko chart bar identification with option to use ATR or traditional method (constant box size)
    if use_atr
        bar_id = ones(Int, size(hlc,1))
        box_sizes = atr(hlc, n=n)
        x = hlc[:,3]
        ref_pt = x[1]
        @inbounds for i in 2:size(x,1)
            if abs(x[i]-ref_pt) >= box_sizes[i]
                ref_pt = x[i]
                bar_id[i:end] .+= 1
            end
        end
        return bar_id
    else
        return renko(hlc[:,3], box_size=box_size)
    end
end

ext: jl


text:
push!(Libdl.DL_LOAD_PATH, abspath(@__FILE__, "..", "..", "deps", "lib"));

if (haskey(ENV, "includelocalwrapper") && ENV["includelocalwrapper"]=="on")
  haskey(ENV, "TestingPath") || (ENV["TestingPath"] = abspath(@__FILE__, "..", "..", "test/CoolProp.jl"));
  if (isfile(ENV["TestingPath"]))
    include(ENV["TestingPath"]);
    warn("** TestingMode is enabeled for: " * ENV["TestingPath"] * " **");
  else
    warn("** $(ENV["TestingPath"]) not exists **");
  end
else
  include("../deps/lib/CoolProp.jl");
end


text:
using Base: @propagate_inbounds

using OffsetArrays

using Oceananigans.Architectures
using Oceananigans.Utils
using Oceananigans.Grids: total_length, interior_indices, interior_parent_indices

import Base: size, length, iterate, getindex, setindex!, lastindex
import Adapt

import Oceananigans.Architectures: architecture
import Oceananigans.Utils: datatuple
import Oceananigans.Grids: total_size, topology, nodes, xnodes, ynodes, znodes, xnode, ynode, znode

@hascuda using CuArrays

"""
    AbstractField{X, Y, Z, A, G}

Abstract supertype for fields located at `(X, Y, Z)` with data stored in a container
of type `A`. The field is defined on a grid `G`.
"""
abstract type AbstractField{X, Y, Z, A, G} end

"""
    Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}

A field defined at the location (`X`, `Y`, `Z`), each of which can be either `Cell`
or `Face`, and with data stored in a container of type `A` (typically an array).
The field is defined on a grid `G` and has field boundary conditions `B`.
"""
struct Field{X, Y, Z, A, G, B} <: AbstractField{X, Y, Z, A, G}
                   data :: A
                   grid :: G
    boundary_conditions :: B

    function Field{X, Y, Z}(data, grid, bcs) where {X, Y, Z}
        Tx, Ty, Tz = total_size((X, Y, Z), grid)

        if size(data) != (Tx, Ty, Tz)
            e = "Cannot construct field at ($X, $Y, $Z) with size(data)=$(size(data)). " *
                "`data` must have size ($Tx, $Ty, $Tz)."
            throw(ArgumentError(e))
        end

        return new{X, Y, Z, typeof(data), typeof(grid), typeof(bcs)}(data, grid, bcs)
    end
end

"""
    Field(X, Y, Z, arch, grid, [  bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
                                 data = zeros(arch, grid, (X, Y, Z)) ] )

Construct a `Field` on `grid` with `data` on architecture `arch` with
boundary conditions `bcs`. Each of `(X, Y, Z)` is either `Cell` or `Face` and determines
the field's location in `(x, y, z)`.

Example
=======

julia> ω = Field(Face, Face, Cell, CPU(), RegularCartesianmodel.grid)

"""
function Field(X, Y, Z, arch, grid,
                bcs = FieldBoundaryConditions(grid, (X, Y, Z)),
               data = zeros(eltype(grid), arch, grid, (X, Y, Z)))

    return Field{X, Y, Z}(data, grid, bcs)
end

#####
##### Convenience constructor for Field that uses a 3-tuple of locations rather than a list of locations:
#####

# Type "destantiation": convert Face() to Face and Cell() to Cell if needed.
destantiate(X) = typeof(X)
destantiate(X::DataType) = X

"""
    Field(L::Tuple, arch, grid, data, bcs)

Construct a `Field` at the location defined by the 3-tuple `L`,
whose elements are `Cell` or `Face`.
"""
Field(L::Tuple, args...) = Field(destantiate.(L)..., args...)

#####
##### Special constructors for tracers and velocity fields
#####

"""
    CellField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
              [  bcs = TracerBoundaryConditions(grid),
                data = zeros(FT, arch, grid, (Cell, Cell, Cell) ] )

Return a `Field{Cell, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function CellField(FT::DataType, arch, grid,
                    bcs = TracerBoundaryConditions(grid),
                   data = zeros(FT, arch, grid, (Cell, Cell, Cell)))

    return Field{Cell, Cell, Cell}(data, grid, bcs)
end

"""
    XFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
               [  bcs = UVelocityBoundaryConditions(grid),
                 data = zeros(FT, arch, grid, (Face, Cell, Cell) ] )

Return a `Field{Face, Cell, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function XFaceField(FT::DataType, arch, grid,
                     bcs = UVelocityBoundaryConditions(grid),
                    data = zeros(FT, arch, grid, (Face, Cell, Cell)))

    return Field{Face, Cell, Cell}(data, grid, bcs)
end

"""
    YFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
               [  bcs = VVelocityBoundaryConditions(grid),
                 data = zeros(FT, arch, grid, (Cell, Face, Cell)) ] )

Return a `Field{Cell, Face, Cell}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function YFaceField(FT::DataType, arch, grid,
                     bcs = VVelocityBoundaryConditions(grid),
                    data = zeros(FT, arch, grid, (Cell, Face, Cell)))

    return Field{Cell, Face, Cell}(data, grid, bcs)
end

"""
    ZFaceField([ FT=eltype(grid) ], arch::AbstractArchitecture, grid,
               [  bcs = WVelocityBoundaryConditions(grid),
                 data = zeros(FT, arch, grid, (Cell, Cell, Face)) ] )

Return a `Field{Cell, Cell, Face}` on architecture `arch` and `grid` containing `data`
with field boundary conditions `bcs`.
"""
function ZFaceField(FT::DataType, arch, grid,
                     bcs = WVelocityBoundaryConditions(grid),
                    data = zeros(FT, arch, grid, (Cell, Cell, Face)))

    return Field{Cell, Cell, Face}(data, grid, bcs)
end

 CellField(arch::AbstractArchitecture, grid, args...) =  CellField(eltype(grid), arch, grid, args...)
XFaceField(arch::AbstractArchitecture, grid, args...) = XFaceField(eltype(grid), arch, grid, args...)
YFaceField(arch::AbstractArchitecture, grid, args...) = YFaceField(eltype(grid), arch, grid, args...)
ZFaceField(arch::AbstractArchitecture, grid, args...) = ZFaceField(eltype(grid), arch, grid, args...)

#####
##### Functions for querying fields
#####

location(a) = nothing

"Returns the location `(X, Y, Z)` of an `AbstractField{X, Y, Z}`."
location(::AbstractField{X, Y, Z}) where {X, Y, Z} = (X, Y, Z) # note no instantiation
location(f, i) = location(f)[i]

"Returns the architecture where the field data `f.data` is stored."
architecture(f::Field) = architecture(f.data)
architecture(o::OffsetArray) = architecture(o.parent)

"Returns the length of a field's `data`."
@inline length(f::Field) = length(f.data)

"Returns the topology of a fields' `grid`."
@inline topology(f, args...) = topology(f.grid, args...)

"""
    size(loc, grid)

Returns the size of a field at `loc` on `grid`.
This is a 3-tuple of integers corresponding to the number of interior nodes
of `f` along `x, y, z`.
"""
@inline size(loc, grid) = (total_length(loc[1], topology(grid, 1), grid.Nx),
                           total_length(loc[2], topology(grid, 2), grid.Ny),
                           total_length(loc[3], topology(grid, 3), grid.Nz))

"""
    size(f::AbstractField{X, Y, Z}) where {X, Y, Z}

Returns the size of an `AbstractField{X, Y, Z}` located at `X, Y, Z`.
This is a 3-tuple of integers corresponding to the number of interior nodes
of `f` along `x, y, z`.
"""
@inline size(f::AbstractField) = size(location(f), f.grid)

"""
    total_size(loc, grid)

Returns the "total" size of a field at `loc` on `grid`.
This is a 3-tuple of integers corresponding to the number of grid points
contained by `f` along `x, y, z`.
"""
@inline total_size(loc, grid) = (total_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                 total_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                 total_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

@inline total_size(f::AbstractField) = total_size(location(f), f.grid)

@propagate_inbounds getindex(f::Field, inds...) = @inbounds getindex(f.data, inds...)
@propagate_inbounds setindex!(f::Field, v, inds...) = @inbounds setindex!(f.data, v, inds...)

@inline lastindex(f::Field) = lastindex(f.data)
@inline lastindex(f::Field, dim) = lastindex(f.data, dim)

"Returns `f.data` for `f::Field` or `f` for `f::AbstractArray."
@inline data(a) = a # fallback
@inline data(f::Field) = f.data

@inline cpudata(a) = data(a)

@hascuda @inline cpudata(f::Field{X, Y, Z, <:CuArray}) where {X, Y, Z} =
    OffsetArray(Array(parent(f)), f.grid, location(f))

# Endpoint for recursive `datatuple` function:
@inline datatuple(obj::AbstractField) = data(obj)

"Returns `f.data.parent` for `f::Field`."
@inline Base.parent(f::Field) = f.data.parent

"Returns a view of `f` that excludes halo points."
@inline interior(f::Field{X, Y, Z}) where {X, Y, Z} = view(f.data, interior_indices(X, topology(f, 1), f.grid.Nx),
                                                                   interior_indices(Y, topology(f, 2), f.grid.Ny),
                                                                   interior_indices(Z, topology(f, 3), f.grid.Nz))

"Returns a reference (not a view) to the interior points of `field.data.parent.`"
@inline interiorparent(f::Field{X, Y, Z}) where {X, Y, Z} =
    @inbounds f.data.parent[interior_parent_indices(X, topology(f, 1), f.grid.Nx, f.grid.Hx),
                            interior_parent_indices(Y, topology(f, 2), f.grid.Ny, f.grid.Hy),
                            interior_parent_indices(Z, topology(f, 3), f.grid.Nz, f.grid.Hz)]

iterate(f::Field, state=1) = iterate(f.data, state)

@inline xnode(i, ψ::Field{X, Y, Z}) where {X, Y, Z} = xnode(X, i, ψ.grid)
@inline ynode(j, ψ::Field{X, Y, Z}) where {X, Y, Z} = ynode(Y, j, ψ.grid)
@inline znode(k, ψ::Field{X, Y, Z}) where {X, Y, Z} = znode(Z, k, ψ.grid)

xnodes(ψ::AbstractField) = xnodes(location(ψ, 1), ψ.grid)
ynodes(ψ::AbstractField) = ynodes(location(ψ, 2), ψ.grid)
znodes(ψ::AbstractField) = znodes(location(ψ, 3), ψ.grid)

nodes(ψ::AbstractField) = (xnodes(ψ), ynodes(ψ), znodes(ψ))

#####
##### Creating offset arrays for field data by dispatching on architecture.
#####

"""
Return a range of indices for a field located at `Cell` centers
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(loc, topo, N, H=0) = 1 - H : N + H

"""
Return a range of indices for a field located at cell `Face`s
`along a grid dimension of length `N` and with halo points `H`.
"""
offset_indices(::Type{Face}, ::Type{Bounded}, N, H=0) = 1 - H : N + H + 1

"""
    OffsetArray(underlying_data, grid, loc)

Returns an `OffsetArray` that maps to `underlying_data` in memory,
with offset indices appropriate for the `data` of a field on
a `grid` of `size(grid)` and located at `loc`.
"""
function OffsetArray(underlying_data, grid, loc)
    ii = offset_indices(loc[1], topology(grid, 1), grid.Nx, grid.Hx)
    jj = offset_indices(loc[2], topology(grid, 2), grid.Ny, grid.Hy)
    kk = offset_indices(loc[3], topology(grid, 3), grid.Nz, grid.Hz)

    return OffsetArray(underlying_data, ii, jj, kk)
end

"""
    zeros([FT=Float64], ::CPU, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT`, with
parent data in CPU memory and indices corresponding to a field on a
`grid` of `size(grid)` and located at `loc`.
"""
function Base.zeros(FT, ::CPU, grid, loc=(Cell, Cell, Cell))
    underlying_data = zeros(FT, total_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                total_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                total_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    return OffsetArray(underlying_data, grid, loc)
end

"""
    zeros([FT=Float64], ::GPU, grid, loc)

Returns an `OffsetArray` of zeros of float type `FT`, with
parent data in GPU memory and indices corresponding to a field on a `grid`
of `size(grid)` and located at `loc`.
"""
function Base.zeros(FT, ::GPU, grid, loc=(Cell, Cell, Cell))
    underlying_data = CuArray{FT}(undef, total_length(loc[1], topology(grid, 1), grid.Nx, grid.Hx),
                                         total_length(loc[2], topology(grid, 2), grid.Ny, grid.Hy),
                                         total_length(loc[3], topology(grid, 3), grid.Nz, grid.Hz))

    underlying_data .= 0 # Ensure data is initially 0.

    return OffsetArray(underlying_data, grid, loc)
end

# Default to type of Grid
Base.zeros(arch, grid, loc=(Cell, Cell, Cell)) = zeros(eltype(grid), arch, grid, loc)

Base.zeros(FT, ::CPU, Nx, Ny, Nz) = zeros(FT, Nx, Ny, Nz)
Base.zeros(FT, ::GPU, Nx, Ny, Nz) = zeros(FT, Nx, Ny, Nz) |> CuArray

Base.zeros(arch, grid, Nx, Ny, Nz) = zeros(eltype(grid), arch, Nx, Ny, Nz)


text:
using Graphs

abstract AbstractSpikingNeuron
abstract AbstractSpikingLayer
abstract AbstractConnection
abstract AbstractNet

type SimpleSpikingNeuron <: AbstractSpikingNeuron
	lastSpike::Float64
    didSpike::Bool
    function SimpleSpikingNeuron()
    	new(-Inf,false)
	end
end

type HHNeuron <: AbstractSpikingNeuron
	m::Float64 # Na channel activation
	h::Float64 # Na channel inactivation
	n::Float64 # K channel activation
	v::Float64 # membrane potential
	lastSpike::Float64 # time of last spike
	didSpike::Bool
end

typealias SynapseType Int64
const EXC = 1
const INH = 2

type SpikingSynapse
	stype::SynapseType
	s::Int64 # source index
	t::Int64 # target index
	w::Float64 # weight
	d::Float64 # delay
	g::Float64 # transmitted conductance
	tau::Float64 # synaptic time constant
	incoming::Vector{Bool}
	_head::Int64
	_ql::Int64
end

type SpikingLayer{T <: AbstractSpikingNeuron} <: AbstractSpikingLayer
	name::ASCIIString
	N::Int64
	units::Vector{T}
end

type SpikingConnection <: AbstractConnection
	s::Int64 # source layer index (pre)
	t::Int64 # target layer index (post)
	synapses::Vector{SpikingSynapse}
	stdp_enabled::Bool
	avg::Float64 # average weight of all synapses
end

type RunSettings
    T # time
    stim_spikes::Vector{Any} # dense matrix or nothing
    # stim_inj::Dict{Int64,Array{Float64,2}}
#    stim::Function  # input(a,b,c) = a'th neuron at timestep b, layer c
end

type RunResults
    names::Dict{Int64, ASCIIString}
    N::Dict{Int64, Int64}
    spikes::Dict{Int64, Vector{(Int64,Int64)}} # sparse spike arrays = time, index
    # pre_init::Vector{Vector{SpikingConnection}} # initial synapses & weights
    # pre_end::Vector{Vector{SpikingConnection}} # final synapses & weights
    dt::Float64
    T::Float64
end

# for building smaller, simple networks
# no concept of layers, just neurons wired together in a small graph.
type SpikingNet <: AbstractNet
	g::SimpleGraph
	l::Vector{AbstractSpikingLayer}
	synapses::Vector{SpikingSynapse}
	stdp_enabled::Bool
end

type SpikingLayersNet <: AbstractNet
    g::SimpleGraph # we don't actually use this graph for anything... except maybe drawing?
    l::Vector{AbstractSpikingLayer} # layers, represented by each node
    # SpikingConnection in pre and post are modified by reference!
    pre::Vector{Vector{SpikingConnection}} # presynaptic synapses for each layer
    post::Vector{Vector{SpikingConnection}} # postsynaptic synapses for each layer
end


text:
module ColonyRuntime
include("./core.jl")
include("./rpc.jl")
include("./Crypto.jl")

import .Crypto

using JSON
using JSON2
using HTTP
using Base64

struct ColoniesServer
  host::String
  port::Int64
end

function addcolony(server::ColoniesServer, colony::Colony, prvkey::String)
  rpcmsg = AddColonyRPC(colony, "addcolonymsg")
  rpcjson = JSON2.write(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "addcolonymsg", payload)

  payload, payloadtype = sendrpcmsg(rpcmsg)
  JSON2.read(payload, Colony)
end

function addruntime(server::ColoniesServer, runtime::Runtime, prvkey::String)
  rpcmsg = AddRuntimeRPC(runtime, "addruntimemsg")
  rpcjson = JSON.json(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "addruntimemsg", payload)

  payload, payloadtype = sendrpcmsg(rpcmsg)
  JSON2.read(payload, Runtime)
end

function approveruntime(server::ColoniesServer, runtimeid::String, prvkey::String)
  rpcmsg = ApproveRuntimeRPC(runtimeid, "approveruntimemsg")
  rpcjson = JSON.json(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "approveruntimemsg", payload)

  sendrpcmsg(rpcmsg)
end

function submitprocess(server::ColoniesServer, spec::ProcessSpec, prvkey::String)
  rpcmsg = SubmitProcessSpecRPC(spec, "submitprocessespecmsg")
  rpcjson = JSON.json(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "submitprocessespecmsg", payload)

  payload, payloadtype = sendrpcmsg(rpcmsg)
  JSON2.read(payload, Process)
end

function getprocess(server::ColoniesServer, processid::String, prvkey::String)
  rpcmsg = GetProcessRPC(processid, "getprocessmsg")
  rpcjson = JSON.json(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "getprocessmsg", payload)

  payload, payloadtype = sendrpcmsg(rpcmsg)
  JSON2.read(payload, Process)
end

function assignprocess(server::ColoniesServer, colonyid::String, prvkey::String)
  rpcmsg = AssignProcessRPC(colonyid, "assignprocessmsg")
  rpcjson = JSON.json(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "assignprocessmsg", payload)

  payload, payloadtype = sendrpcmsg(rpcmsg)
  JSON2.read(payload, Process)
end

function addattribute(server::ColoniesServer, attribute::Attribute, prvkey::String)
  rpcmsg = AddAttributeRPC(attribute, "addattributemsg")
  rpcjson = JSON.json(rpcmsg)

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, "addattributemsg", payload)

  payload, payloadtype = sendrpcmsg(rpcmsg)
  JSON2.read(payload, Attribute)
end

function closeprocess(server::ColoniesServer, processid::String, successful::Bool, prvkey::String)
  payloadtype = "closesuccessfulmsg"
  rpcmsg = CloseSuccessfulRPC(processid, payloadtype)
  rpcjson = JSON.json(rpcmsg)

  if !successful
    payloadtype = "closefailedmsg"
    rpcmsg = CloseSuccessfulRPC(processid, payloadtype)
    rpcjson = JSON.json(rpcmsg)
  end

  payload = base64enc(rpcjson)
  sig = Crypto.sign(payload, prvkey)
  rpcmsg = RPCMsg(sig, payloadtype, payload)

  sendrpcmsg(rpcmsg)
end

end


text:
using ITensors

let
  Ny = 6
  Nx = 12

  N = Nx*Ny

  sites = siteinds("S=1/2",N;conserve_qns=true)

  # Turning on QN conservation can
  # give large speedups for DMRG:
  #sites = siteinds("S=1/2",N,conserve_qns=true)

  lattice = square_lattice(Nx,Ny,yperiodic=false)

  ampo = AutoMPO()
  for b in lattice
    ampo += (0.5,"S+",b.s1,"S-",b.s2)
    ampo += (0.5,"S-",b.s1,"S+",b.s2)
    ampo += ("Sz",b.s1,"Sz",b.s2)
  end
  H = MPO(ampo,sites)

  state = [isodd(n) ? "Up" : "Dn" for n=1:N]
  # Initialize wavefunction to a random MPS
  # of bond-dimension 10 with same quantum
  # numbers as `state`
  psi0 = randomMPS(sites,state,20)

  sweeps = Sweeps(10)
  maxdim!(sweeps,20,60,100,100,200,400,800)
  cutoff!(sweeps,1E-8)
  @show sweeps

  energy,psi = dmrg(H,psi0,sweeps)

  return
end


text:
import VimBindings.Motion

@testset "motion addition" begin
    m1 = Motion(0, 4)
    m2 = Motion(4, 8)

    @test m1 + m2 == Motion(0, 8)

    @test Motion(1, 3) + Motion(3, 5) == Motion(1, 5)

    @test Motion(3, 0) + Motion(3, 5) == Motion(0, 5)
end


text:
using Base.Test
using NGSIM
using NBInclude

nbinclude(joinpath(dirname(@__FILE__), "..", "jnotebooks", "Demo.ipynb"))


text:
push!(LOAD_PATH, ".")
using WildBootTests
using StatFiles, StatsModels, DataFrames, DataFramesMeta, BenchmarkTools, Plots, CategoricalArrays, Random, StableRNGs

open("unittests.log", "w") do log  # use Github Desktop to detect changes in output

df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\collapsed.dta"))
dropmissing!(df)
f = @formula(hasinsurance ~ 1 + selfemployed + post + post_self)
f = apply_schema(f, schema(f, df, Dict(:hasinsurance => ContinuousTerm)))
resp, predexog = modelcols(f, df)

println(log, "\nregress hasinsurance selfemployed post post_self, cluster(year)")
println(log, "\nboottest post_self=.04, weight(webb)")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), auxwttype=:webb, rng=StableRNG(1231))
println(log, test)

test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), auxwttype=:webb, rng=StableRNG(1231), issorted=true)
println(log, test)

println(log, "\nboottest post_self=.04, weight(webb) reps(9999999) noci")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), reps=9999999, auxwttype=:webb, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nboottest post_self=.04, weight(normal) reps(9999) noci")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:normal, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nboottest post_self=.04, weight(gamma) reps(9999) noci svv")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:gamma, getci=false, rng=StableRNG(1231), getauxweights=true)
println(log, test)

println(log, "\nboottest post_self=.04, weight(mammen) reps(9999) noci")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:mammen, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nboottest post_self=.04, weight(mammen) reps(9999) boottype(score)")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:mammen, scorebs=true, rng=StableRNG(1231))
println(log, test)

println(log, "\nregress hasinsurance selfemployed post post_self, robust")
println(log, "\nboottest post_self=.04, weight(webb)")
test = wildboottest([0 0 0 1], [.04]; resp, predexog, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)

println(log, "\nregress hasinsurance selfemployed post post_self, cluster(year)")
println(log, "boottest (post_self=.05) (post=-.02), reps(9999) weight(webb)")
test = wildboottest([0 0 0 1; 0 0 1 0], [.05; -.02]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)

test = wildboottest([0 0 0 1; 0 0 1 0], [.05; -.02]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:webb, rng=StableRNG(1231), bootstrapc=true)
println(log, test)

println(log, "boottest (post_self=.05) (post=-.02) (selfemployed=-.15), reps(9999) weight(webb)")
test = wildboottest([0 0 0 1; 0 0 1 0; 0 1 0 0], [.05; -.02; -.15]; resp, predexog, clustid=Int32.(df.year), reps=9999, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)

println(log, "\nregress hasinsurance selfemployed post post_self")
println(log, "\nboottest post_self=.04, weight(webb)")
test = wildboottest([0 0 0 1], [.04]               ; resp, predexog, hetrobust=false, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)
println(log, "boottest (post_self=.05) (post=-.02), reps(999) weight(webb)")  # bad
test = wildboottest([0 0 0 1; 0 0 1 0], [.05; -.02]; resp, predexog, hetrobust=false, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)
println(log, "scoretest (post_self=.05)")
test = wildboottest([0 0 0 1], [.05]; resp, predexog, hetrobust=false, scorebs=true, reps=0)
println(log, test)
println(log, "scoretest (post_self=.05) (post=-.02)")
test = wildboottest([0 0 0 1; 0 0 1 0], [.05; -.02]; resp, predexog, hetrobust=false, scorebs=true, reps=0)
println(log, test)
println(log, "boottest (post_self=.08), boottype(score)")
test = wildboottest([0 0 0 1], [.08]; resp, predexog, hetrobust=false, scorebs=true, rng=StableRNG(1231))
println(log, test)
println(log, "boottest (post_self=.05) (post=-.02), boottype(score)")
test = wildboottest([0 0 0 1; 0 0 1 0], [.05; -.02]; resp, predexog, hetrobust=false, scorebs=true, rng=StableRNG(1231))
println(log, test)

df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))[:,[:wage; :tenure; :ttl_exp; :collgrad; :industry]]
dropmissing!(df)
desc = describe(df, :eltype)
f = @formula(wage ~ 1 + tenure + ttl_exp + collgrad)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

println(log, "\nconstraint 1 ttl_exp = .2")
println(log, "cnsreg wage tenure ttl_exp collgrad, constr(1) cluster(industry)")
println(log, "boottest tenure")
test = wildboottest([0 1 0 0], [0]; R1=[0 0 1 0], r1=[.2], resp, predexog, clustid=df.industry, rng=StableRNG(1231))
println(log, test)

test = wildboottest([0 1 0 0], [0]; R1=[0 0 1 0], r1=[.2], resp, predexog, clustid=[collect(1:1000); collect(1:1217)], reps=9999, rng=StableRNG(1231))  # granular but not pure robust
println(log, test)

test = wildboottest([0 1 0 0], [0]; resp, predexog, clustid=[collect(1:1000); collect(1:1217)], feid=df.industry, reps=9999, rng=StableRNG(1231))
println(log, test)

println(log, "\nivregress liml wage ttl_exp collgrad (tenure = union), cluster(industry)")
println(log, "boottest tenure, ptype(equaltail) reps(9999)")

df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage; :tenure; :ttl_exp; :collgrad; :industry; :union]]
dropmissing!(df)
f = @formula(wage ~ 1 + ttl_exp + collgrad)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
ivf = @formula(tenure ~ union)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=9999, ptype=:equaltail, rng=StableRNG(1231))
println(log, test)

println(log, "boottest tenure, nonull reps(9999) matsize(.1)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=99999, imposenull=false, rng=StableRNG(1231), maxmatsize=.1)
println(log, test)

println(log, "boottest tenure, ptype(upper) svmat(t) reps(9999)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=9999, ptype=:upper, diststat=:t, rng=StableRNG(1231))
println(log, test)

println(log, "boottest tenure, ptype(lower) svmat(numer) reps(9999)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=9999, ptype=:lower, diststat=:numer, rng=StableRNG(1231))
println(log, test)

println(log, "\nivregress liml wage ttl_exp collgrad (tenure = union), cluster(industry)")
println(log, "boottest tenure, ptype(equaltail) reps(9999)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, liml=true, reps=9999, ptype=:equaltail, rng=StableRNG(1231))
println(log, test)

println(log, "\nivregress liml wage ttl_exp collgrad (tenure = union) if industry<., robust")
println(log, "boottest tenure, ptype(equaltail) reps(99) noci")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, small=false, liml=true, reps=99, ptype=:equaltail, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nivregress 2sls wage ttl_exp collgrad (tenure = union) if industry<., robust")
println(log, "boottest tenure, ptype(equaltail) reps(99) noci")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, small=false, reps=99, ptype=:equaltail, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "boottest collgrad tenure, ptype(equaltail) reps(99) noci")
test = wildboottest([0 0 0 1; 0 0 1 0], [0;0]; resp, predexog, predendog, inst, small=false, reps=99, ptype=:equaltail, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nivregress 2sls wage ttl_exp collgrad (tenure = union) if industry<.")
println(log, "boottest tenure, ptype(equaltail) reps(99) noci")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, hetrobust=false, small=false, reps=99, ptype=:equaltail, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "boottest tenure collgrad, ptype(equaltail) reps(99) noci")
test = wildboottest([0 0 0 1; 0 0 1 0], [0;0]; resp, predexog, predendog, inst, hetrobust=false, small=false, reps=99, ptype=:equaltail, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nivregress 2sls wage ttl_exp collgrad (tenure = union), cluster(industry)")
println(log, "\nboottest tenure, ptype(equaltail) weight(webb) stat(c) gridmin(-5) gridmax(5) gridpoints(100)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, auxwttype=:webb, bootstrapc=true, ptype=:equaltail, rng=StableRNG(1231), gridmin=[-5], gridmax=[5], gridpoints=[100])
println(log, test)

println(log, "\nboottest tenure, ptype(equaltail) weight(webb) stat(c) gridmin(-5) gridmax(5) gridpoints(100) maxmatsize(.01)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, auxwttype=:webb, bootstrapc=true, ptype=:equaltail, rng=StableRNG(1231), gridmin=[-5], gridmax=[5], gridpoints=[100], maxmatsize=.01)
println(log, test)

println(log, "\nivregress 2sls wage ttl_exp collgrad (tenure = union) if industry<., robust")
println(log, "\nboottest tenure, ptype(equaltail) maxmatsize(0.005) noci weight(webb) ptype(equaltail)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, small=false, auxwttype=:webb, ptype=:equaltail, rng=StableRNG(1231), getci=false, maxmatsize=.005)
println(log, test)

test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=[collect(1:1000); collect(1:855)], small=false, auxwttype=:webb, ptype=:equaltail, getci=false, rng=StableRNG(1231), maxmatsize=.005)
println(log, test)

println(log, "\nivregress 2sls wage ttl_exp collgrad (tenure = union), cluster(industry)")
println(log, "\nboottest, ar")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, arubin=true, rng=StableRNG(1231))
println(log, test)
plot(plotpoints(test)...)

println(log, "\nboottest, ar nonull")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, arubin=true, rng=StableRNG(1231), imposenull=false)
println(log, test)

println(log, "\nscoretest tenure")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=0)
println(log, test)

println(log, "\nwaldtest tenure, pytpe(upper)")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, small=false, reps=0, imposenull=false, getplot=false, ptype=:upper)
println(log, test)

println(log, "\nivregress liml wage (tenure = collgrad ttl_exp), cluster(industry)")
println(log, "boottest tenure")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage, :tenure, :ttl_exp, :collgrad, :industry]]
dropmissing!(df)
f = @formula(wage ~ 1)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
ivf = @formula(tenure ~ collgrad + ttl_exp)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
test = wildboottest([0 1], [0]; resp, predexog, predendog, inst, liml=true, clustid=df.industry, small=false, rng=StableRNG(1231))
println(log, test)

println(log, "\nivreg2 wage collgrad smsa race age (tenure = union married), cluster(industry) fuller(1)")
println(log, "boottest tenure, nograph weight(webb) reps(9999)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
df = df[:, [:wage, :tenure, :ttl_exp, :collgrad, :smsa, :race, :age, :union, :married, :industry]]
dropmissing!(df)
f = @formula(wage ~ 1 + collgrad + smsa + race + age)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
ivf = @formula(tenure ~ union + married)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
test = wildboottest([0 0 0 0 0 1], [0]; resp, predexog, predendog, inst, fuller=1, clustid=df.industry, small=false, reps=9999, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)

println(log, "boottest tenure, noci bootcluster(individual) weight(webb)")
test = wildboottest([0 0 0 0 0 1], [0]; resp, predexog, predendog, inst, fuller=1, clustid=[collect(1:nrow(df)) df.industry], nbootclustvar=1, nerrclustvar=1, small=false, auxwttype=:webb, getci=false, rng=StableRNG(1231))
println(log, test)

println(log, "boottest tenure, nograph bootcluster(collgrad) cluster(collgrad industry) weight(webb) reps(9999)")
test = wildboottest([0 0 0 0 0 1], [0]; resp, predexog, predendog, inst, fuller=1, clustid=clustid=Matrix(df[:, [:collgrad, :industry]]), nbootclustvar=1, nerrclustvar=2, small=false, reps=9999, auxwttype=:webb, rng=StableRNG(1231))
println(log, test)

println(log, "\nareg wage ttl_exp collgrad tenure [aw=hours] if occupation<., cluster(age) absorb(industry)")
println(log, "boottest tenure, cluster(age occupation) bootcluster(occupation)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\nlsw88.dta"))
dropmissing!(df)
f = @formula(wage ~ ttl_exp + collgrad + tenure)  # constant unneeded in FE model
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
test = wildboottest([0 0 1], [0]; resp, predexog, clustid=Matrix(df[:, [:occupation, :age]]), nbootclustvar=1, nerrclustvar=2, obswt=df.hours, feid=df.industry, rng=StableRNG(1231))
println(log, test)

println(log, "\nareg wage ttl_exp collgrad tenure if occupation<., robust absorb(industry)")
println(log, "boottest tenure")
test = wildboottest([0 0 1], [0]; resp, predexog, feid=df.industry, rng=StableRNG(1231))
println(log, test)

println(log, "\nareg wage ttl_exp collgrad tenure [aw=hours] if occupation<., robust absorb(industry)")
println(log, "boottest tenure")
test = wildboottest([0 0 1], [0]; resp, predexog, obswt=df.hours, feid=df.industry, rng=StableRNG(1231))
println(log, test)

println(log, "\nivreghdfe wage ttl_exp collgrad tenure (occupation = union married) [aw=hours], liml cluster(industry) absorb(industry)")
ivf = @formula(occupation ~ union + married)
ivf = apply_schema(ivf, schema(ivf, df))
predendog, inst = modelcols(ivf, df)
println(log, "boottest tenure")
test = wildboottest([0 0 1 0], [0]; resp, predexog, predendog, inst, clustid=df.industry, obswt=df.hours, feid=df.industry, fedfadj=1, rng=StableRNG(1231), liml=true)
println(log, test)
println(log, "boottest occupation")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, obswt=df.hours, feid=df.industry, fedfadj=1, rng=StableRNG(1231), liml=true)
println(log, test)

println(log, "\nivreghdfe wage ttl_exp collgrad tenure (occupation = union married) [aw=hours] if grade<., liml cluster(industry) absorb(age)")
println(log, "boottest tenure")
test = wildboottest([0 0 1 0], [0]; resp, predexog, predendog, inst, clustid=df.industry, obswt=df.hours, feid=df.age, rng=StableRNG(1231), liml=true)
println(log, test)
println(log, "boottest collgrad tenure")
test = wildboottest([0 0 1 0; 0 1 0 0], [0; 0]; resp, predexog, predendog, inst, clustid=df.industry, obswt=df.hours, feid=df.age, rng=StableRNG(1231), liml=true, reps=99)
println(log, test)
println(log, "boottest occupation")
test = wildboottest([0 0 0 1], [0]; resp, predexog, predendog, inst, clustid=df.industry, obswt=df.hours, feid=df.age, rng=StableRNG(1231), liml=true, gridmin=[-1], gridmax=[1])
println(log, test)

println(log, "boottest tenure | _b[collgrad] = 0")
test = wildboottest([0 0 1 0], [0]; R1=[0 1 0 0], r1=[0], resp, predexog                  , predendog, inst                     , clustid=df.industry, obswt=df.hours, feid=df.age, rng=StableRNG(1231), liml=true, gridmin=[-1], gridmax=[1])
println(log, test)
test = wildboottest([0 1 0], [0];                         resp, predexog=predexog[:,[1,3]], predendog, inst=[inst predexog[:,2]], clustid=df.industry, obswt=df.hours, feid=df.age, rng=StableRNG(1231), liml=true, gridmin=[-1], gridmax=[1])
println(log, test)

df = DataFrame(load(raw"d:\OneDrive\Documents\Macros\abdata.dta"))[:,[:n; :w; :k; :ys; :id; :year; :ind]]
dropmissing!(df)
f = @formula(n ~ w + k)  # constant unneeded in FE model
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)
println(log, "\nareg n w k, absorb(ind)")
println(log, "\nboottest k, cluster(id year)")
test = WildBootTests.wildboottest([0 1], [0]; resp, predexog, clustid=round.(Int,Matrix(df[:, [:id, :year]])), nbootclustvar=2, nerrclustvar=2, feid=round.(Int,df.ind), rng=StableRNG(1231))
println(log, test)
println(log, "\nareg n w k [aw=ys], absorb(ind)")
println(log, "\nboottest k, cluster(id year)")
test = WildBootTests.wildboottest([0 1], [0]; resp, predexog, clustid=round.(Int,Matrix(df[:, [:id, :year]])), nbootclustvar=2, nerrclustvar=2, feid=round.(Int,df.ind), rng=StableRNG(1231), obswt=df.ys)
println(log, test)

println(log, "\nglobal pix lnkm pixpetro pixdia pixwaterd pixcapdist pixmal pixsead pixsuit pixelev pixbdist")
println(log, "global geo lnwaterkm lnkm2split mean_elev mean_suit malariasuit petroleum diamondd")
println(log, "global poly capdistance1 seadist1 borderdist1")
println(log, "encode pixwbcode, gen(ccode)  // make numerical country identifier")
println(log, "qui areg lnl0708s centr_tribe lnpd0 \$pix \$geo \$poly, absorb(ccode)")
println(log, "boottest centr_tribe, nogr reps(9999) clust(ccode pixcluster) bootcluster(ccode)")
println(log, "boottest centr_tribe, nogr reps(9999) clust(ccode pixcluster) bootcluster(pixcluster)")
println(log, "boottest centr_tribe, nogr reps(9999) clust(ccode pixcluster) bootcluster(ccode pixcluster)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Work\Econometrics\Wild cluster\pixel-level-baseline-final.dta"))
pix  = [:lnkm, :pixpetro, :pixdia, :pixwaterd, :pixcapdist, :pixmal, :pixsead, :pixsuit, :pixelev, :pixbdist]
geo  = [:lnwaterkm, :lnkm2split, :mean_elev, :mean_suit, :malariasuit, :petroleum, :diamondd]
poly = [:capdistance1, :seadist1, :borderdist1]
df = df[:,[pix; geo; poly; :lnl0708s; :centr_tribe; :lnpd0; :pixwbcode; :pixcluster]]
dropmissing!(df)
df.ccode = levelcode.(categorical(df.pixwbcode, compress=true))
df.pixcode = levelcode.(categorical(df.pixcluster, compress=true))
f = Term(:lnl0708s) ~ sum(term.([:centr_tribe; :lnpd0; pix; geo; poly]))
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

test = wildboottest([1 zeros(1,size(predexog,2)-1)], [0]; resp, predexog, clustid=Matrix(df[:, [:ccode, :pixcode]]), nbootclustvar=1, nerrclustvar=2, feid=df.ccode, reps=9999, rng=StableRNG(1231))
println(log, test)
test = wildboottest([1 zeros(1,size(predexog,2)-1)], [0]; resp, predexog, clustid=Matrix(df[:, [:pixcode, :ccode]]), nbootclustvar=1, nerrclustvar=2, feid=df.ccode, reps=9999, rng=StableRNG(1231))
println(log, test)
test = wildboottest([1 zeros(1,size(predexog,2)-1)], [0]; resp, predexog, clustid=Matrix(df[:, [:pixcode, :ccode]]), nbootclustvar=2, nerrclustvar=2, feid=df.ccode, reps=9999, rng=StableRNG(1231))
println(log, test)

println(log, "\ninfile coll merit male black asian year state chst using regm.raw, clear")
println(log, "qui regress coll merit male black asian i.year i.state if !inlist(state,34,57,59,61,64,71,72,85,88), cluster(state)	")
println(log, "generate individual = _n  // unique ID for each observation")
println(log, "boottest merit, nogr reps(9999) gridpoints(10)  // defaults to bootcluster(state)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) bootcluster(state year)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull bootcluster(state year)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) bootcluster(individual)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull bootcluster(individual)")
println(log, "boottest merit, nogr reps(9999) gridpoints(10) nonull bootcluster(individual) maxmatsize(.1)")
df = DataFrame(load(raw"d:\OneDrive\Documents\Work\Econometrics\Wild cluster\regm.dta"))
df = DataFrame(coll=Bool.(df.coll), merit=Bool.(df.merit), male=Bool.(df.male), black=Bool.(df.black), asian=Bool.(df.asian), state=categorical(Int8.(df.state)), year=categorical(Int16.(df.year)))
dropmissing!(df)
df = df[df.state .∉ Ref([34,57,59,61,64,71,72,85,88]),:]
f = @formula(coll ~ 1 + merit + male + black + asian + year + state)
f = apply_schema(f, schema(f, df))
resp, predexog = modelcols(f, df)

test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=levelcode.(df.state), gridpoints=[10], reps=9999, rng=StableRNG(1231))
println(log, test)
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=levelcode.(df.state), reps=9999, imposenull=false, rng=StableRNG(1231))
println(log, test)
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=[levelcode.(df.year) levelcode.(df.state)], nbootclustvar=2, nerrclustvar=1, reps=9999, rng=StableRNG(1231))
println(log, test)
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid=[levelcode.(df.year) levelcode.(df.state)], nbootclustvar=2, nerrclustvar=1, reps=9999, imposenull=false, rng=StableRNG(1231))
println(log, test)
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid= [collect(1:nrow(df)) levelcode.(df.state)], nbootclustvar=1, nerrclustvar=1, reps=9999, rng=StableRNG(1231))
println(log, test)
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid= [collect(1:nrow(df)) levelcode.(df.state)], nbootclustvar=1, nerrclustvar=1, reps=9999, imposenull=false, rng=StableRNG(1231))
println(log, test)
test = wildboottest([0 1 zeros(1,size(predexog,2)-2)], [0]; resp, predexog, clustid= [collect(1:nrow(df)) levelcode.(df.state)], nbootclustvar=1, nerrclustvar=1, reps=9999, imposenull=false, maxmatsize=.1, rng=StableRNG(1231))
println(log, test)
end


text:
struct RO_raw_data{T<:CSV.File}
    Weather_RO::T
    GPP_data_RO::T
    SWC_data_RO::T
end

function Load_RO_data(;weather_file::String="Weather_RO")
    Weather_RO = CSV.File("./Data/"*weather_file*".csv")
    GPP_data_RO = CSV.File("./Data/GPP_RO.csv")
    SWC_data_RO = CSV.File("./Data/Daily_SWC_RO.csv")

    RO_data = RO_raw_data(Weather_RO,GPP_data_RO,SWC_data_RO)

    return RO_data
end

function create_acclimation(;X0 = -4.0,τ =  7.0,Smax = 16.0,Weather_RO_Temp)
    #Weather_RO = CSV.File("./Data/Weather_RO.csv")
    #Weather_RO_Temp = Weather_RO[8767:end]

    #photosynthesis that accounts for temperature acclimation (Mäkelä et al., 2004; Mäkelä et al., 2008).
    #X0 = -4.0 #minium tempreture for Photosynthesis
    #τ =  7.0 #Days
    #Smax = 16.0
    X_vec = Float64[]
    X_growth_vec = Float64[]
    X_growth_week_vec = Float64[]
    Dates_growth_vec = Dates.DateTime[]

    push!(X_vec,Weather_RO_Temp[1][2])
    if Weather_RO_Temp[1][11]==1
        push!(X_growth_vec,Weather_RO_Temp[1][2])
        push!(Dates_growth_vec,Weather_RO_Temp[1][1])
    end
    for i = 2:length(Weather_RO_Temp)
        X = X_vec[i-1]+(Weather_RO_Temp[i][2]-X_vec[i-1])/τ
        push!(X_vec,X)
        if Weather_RO_Temp[i][11]==1
            push!(X_growth_vec,X)
            push!(Dates_growth_vec,Weather_RO_Temp[i][1])
        end
    end

    S_growth_vec = max.(X_growth_vec.-X0,Ref(0))
    fs_growth_vec = min.(S_growth_vec/Smax,Ref(1))
    S_growth_week_vec = Float64[]
    fs_growth_week_vec = Float64[]

    i = 1
    while i <=length(X_growth_vec)-6
        if Dates.year(Dates_growth_vec[i])==Dates.year(Dates_growth_vec[i+6])
            mean_val = mean(X_growth_vec[i:i+6])
            push!(X_growth_week_vec,mean_val)
            push!(S_growth_week_vec,mean(S_growth_vec[i:i+6]))
            push!(fs_growth_week_vec,mean(fs_growth_vec[i:i+6]))
            i += 7
        else
            i += 1
        end
    end

    return fs_growth_week_vec
end

function create_weather_struct_RO(RO_data::RO_raw_data;stand_type::String="Fertilized",X0 = -4.0,τ =  7.0,Smax = 16.0,
    weather_file::String="Weather_RO")
    #Creates weekly mean time series for Rosinedal 2015-2018 (Or 2015-2019 if weather_file=Weather_RO2)
    if stand_type == "Fertilized"
        GPP_ind = 2
        SWC_data_RO_ind = 3
    elseif stand_type == "Control"
        GPP_ind = 3
        SWC_data_RO_ind = 4
    else
        error("Wrong input. Either \"Fertilized\" or \"Control\"")
    end

    Weather_RO_Temp = RO_data.Weather_RO[8767:end]
    SWC_data_RO_Temp =  RO_data.SWC_data_RO[366:end]

    Radiation_vec = Float64[]
    PAR_vec = Float64[]
    VP_vec = Float64[]
    Date_vec = Dates.DateTime[]
    GPP_vec = Float64[]
    Tmean_vec = Float64[]
    θₛ_vec = Float64[]
    daylight_vec = Float64[]

    for i in 1:length(Weather_RO_Temp)
        row = Weather_RO_Temp[i]
        if row[11]==1
            push!(Date_vec,row[1])
            push!(Tmean_vec,row[2])
            push!(GPP_vec,RO_data.GPP_data_RO[i][GPP_ind])
            push!(θₛ_vec,SWC_data_RO_Temp[i][SWC_data_RO_ind]/100)

            val = tryparse(Float64,row[7])
            if val !== nothing
                push!(Radiation_vec,val)
                d = Dates.dayofyear(row[1])
                h = CCPH.daylighthour(64*pi/180,d)*3600
                push!(PAR_vec,val*2.6/h)
                push!(daylight_vec,h)
            else
                println("Radiation data is missing at date: $(row[1]), approximating using  linear regression model")
                β = [0.0003977611889578929; -10.463723925716433]
                d = Dates.dayofyear(row[1])
                h = CCPH.daylighthour(64*pi/180,d)*3600
                Radiation_hat =  β[1]*h+β[2]
                push!(Radiation_vec,Radiation_hat)
                push!(PAR_vec,Radiation_hat*2.3/h)
                push!(daylight_vec,h)
            end

            val = row[8]
            push!(VP_vec,val*100)
            if isnan(val)||ismissing(val)
                println("Vapor deficit data is missing (either NaN oc Missing) at date: $(row[1])")
            end
        end
    end

    VPD_vec = CCPH.SVPₜ.(Tmean_vec).-VP_vec

    Date_week_vec = Dates.DateTime[]
    GPP_week_vec = Float64[]
    PAR_week_vec = Float64[]
    Tmean_week_vec = Float64[]
    VPD_week_vec = Float64[]
    θₛ_week_vec = Float64[]
    daylight_vec = Float64[]

    i = 1
    while i<=length(Date_vec)-6
        if Dates.year(Date_vec[i])==Dates.year(Date_vec[i+6])
            push!(Date_week_vec,Date_vec[i])
            push!(Tmean_week_vec,mean(Tmean_vec[i:i+6]))
            push!(PAR_week_vec,mean(PAR_vec[i:i+6]))
            push!(VPD_week_vec,mean(VPD_vec[i:i+6]))
            push!(θₛ_week_vec,mean(θₛ_vec[i:i+6]))
            push!(GPP_week_vec,mean(GPP_vec[i:i+6]))
            d_vec = Dates.dayofyear.(Date_vec[i:i+6])
            push!(daylight_vec,sum(CCPH.daylighthour.(Ref(64*pi/180),d_vec))*3600)
            i += 7
        else
            i += 1
        end
    end

    acclimation_fac = create_acclimation(;X0 = X0,τ =  τ,Smax = Smax,Weather_RO_Temp)

    # Find the end of each growth period
    if weather_file=="Weather_RO_2"
        t_end_vec = [0,0,0,0,0]
        for j in 1:length(Date_week_vec)
            if Dates.year(Date_week_vec[j]) == 2015
                t_end_vec[1] += 1
                t_end_vec[2] += 1
                t_end_vec[3] += 1
                t_end_vec[4] += 1
                t_end_vec[5] += 1
            elseif Dates.year(Date_week_vec[j]) == 2016
                t_end_vec[2] += 1
                t_end_vec[3] += 1
                t_end_vec[4] += 1
                t_end_vec[5] += 1
            elseif Dates.year(Date_week_vec[j]) == 2017
                t_end_vec[3] += 1
                t_end_vec[4] += 1
                t_end_vec[5] += 1
            elseif Dates.year(Date_week_vec[j]) == 2018
                t_end_vec[4] += 1
                t_end_vec[5] += 1
            else
                t_end_vec[5] += 1
            end
        end

        # calc total daylight for each year
        tot_daylight_hour = [0.0,0.0,0.0,0.0,0.0]
        tot_daylight_hour[1] = sum(daylight_vec[1:t_end_vec[1]])
        tot_daylight_hour[2] = sum(daylight_vec[t_end_vec[1] + 1:t_end_vec[2]])
        tot_daylight_hour[3] = sum(daylight_vec[t_end_vec[2] + 1:t_end_vec[3]])
        tot_daylight_hour[4] = sum(daylight_vec[t_end_vec[3] + 1:t_end_vec[4]])
        tot_daylight_hour[5] = sum(daylight_vec[t_end_vec[4] + 1:t_end_vec[5]])

        # Create total daylight hour time series
        tot_daylight_hour_ts = zeros(length(daylight_vec))
        tot_daylight_hour_ts[1:t_end_vec[1]] .= tot_daylight_hour[1]
        tot_daylight_hour_ts[t_end_vec[1] + 1:t_end_vec[2]] .= tot_daylight_hour[2]
        tot_daylight_hour_ts[t_end_vec[2] + 1:t_end_vec[3]] .= tot_daylight_hour[3]
        tot_daylight_hour_ts[t_end_vec[3] + 1:t_end_vec[4]] .= tot_daylight_hour[4]
        tot_daylight_hour_ts[t_end_vec[4] + 1:t_end_vec[5]] .= tot_daylight_hour[5]
    else
        t_end_vec = [0,0,0,0]
        for j in 1:length(Date_week_vec)
            if Dates.year(Date_week_vec[j]) == 2015
                t_end_vec[1] += 1
                t_end_vec[2] += 1
                t_end_vec[3] += 1
                t_end_vec[4] += 1
            elseif Dates.year(Date_week_vec[j]) == 2016
                t_end_vec[2] += 1
                t_end_vec[3] += 1
                t_end_vec[4] += 1
            elseif Dates.year(Date_week_vec[j]) == 2017
                t_end_vec[3] += 1
                t_end_vec[4] += 1
            else
                t_end_vec[4] += 1
            end
        end

        # calc total daylight for each year
        tot_daylight_hour = [0.0,0.0,0.0,0.0]
        tot_daylight_hour[1] = sum(daylight_vec[1:t_end_vec[1]])
        tot_daylight_hour[2] = sum(daylight_vec[t_end_vec[1] + 1:t_end_vec[2]])
        tot_daylight_hour[3] = sum(daylight_vec[t_end_vec[2] + 1:t_end_vec[3]])
        tot_daylight_hour[4] = sum(daylight_vec[t_end_vec[3] + 1:t_end_vec[4]])

        # Create total daylight hour time series
        tot_daylight_hour_ts = zeros(length(daylight_vec))
        tot_daylight_hour_ts[1:t_end_vec[1]] .= tot_daylight_hour[1]
        tot_daylight_hour_ts[t_end_vec[1] + 1:t_end_vec[2]] .= tot_daylight_hour[2]
        tot_daylight_hour_ts[t_end_vec[2] + 1:t_end_vec[3]] .= tot_daylight_hour[3]
        tot_daylight_hour_ts[t_end_vec[3] + 1:t_end_vec[4]] .= tot_daylight_hour[4]
    end

    weatherts = WeatherTS(Date_week_vec,daylight_vec,tot_daylight_hour_ts,PAR_week_vec,
    Tmean_week_vec,VPD_week_vec,θₛ_week_vec,acclimation_fac)

    return weatherts, GPP_week_vec
end

function plot_weather_struct_RO(weatherts::WeatherTS,GPP_week_vec::Array{Float64,1})
    pl1 = plot(weatherts.date,weatherts.daylight/3600,xlabel="Date",ylabel="Daylight (h)",seriestype=:scatter)
    pl2 = plot(weatherts.date,weatherts.PAR*10^6,xlabel="Date",ylabel="PAR (μmol s⁻¹ m⁻²)",seriestype=:scatter)
    pl3 = plot(weatherts.date,weatherts.temp,xlabel="Date",ylabel="Temp (ᵒC)",seriestype=:scatter)
    pl4 = plot(weatherts.date,weatherts.VPD/1000,xlabel="Date",ylabel="VPD (kPa)",seriestype=:scatter)
    pl5 = plot(weatherts.date,weatherts.θₛ*100,xlabel="Date",ylabel="SWC (%)",seriestype=:scatter)
    pl6 = plot(weatherts.date,CCPH.θₛ2ψₛ.(weatherts.θₛ),xlabel="Date",ylabel="SWP (MPa)",seriestype=:scatter)
    pl7 = plot(weatherts.date,GPP_week_vec,xlabel="Date",ylabel="GPP (g C m⁻² d⁻¹)",seriestype=:scatter)
    pl8 = plot(weatherts.date,weatherts.acclimation_fac,xlabel="Date",ylabel="Acclimation",seriestype=:scatter)

    plot(pl1,pl2,pl3,pl4,pl5,pl6,pl7,pl8,layout=8,legends=false)
end


text:
function sample_mh_jump(model, samples=1000 ; δ=180, σ=5, p=0.2)
    mh = MH(:startpos => x -> [MixtureModel(Normal.([v -  δ, v + δ], σ), [p/2, 1-p, p/2]) for v in x])
    sample(model, mh, samples)
end


text:
using KitBase, FluxReconstruction, OrdinaryDiffEq, LinearAlgebra, Plots
using KitBase.ProgressMeter: @showprogress
using Base.Threads: @threads

begin
    x0 = 0
    x1 = 1
    ncell = 100
    dx = (x1 - x0) / ncell
    deg = 7 # polynomial degree
    nsp = deg + 1
    γ = 5 / 3
    cfl = 0.01
    dt = cfl * dx
    t = 0.0
end
ps = FRPSpace1D(x0, x1, ncell, deg)
ℓ = FR.basis_norm(ps.deg)

u = zeros(ncell, nsp, 3)
for i = 1:ncell, ppp1 = 1:nsp
    if ps.x[i] <= 0.5
        prim = [1.0, 0.0, 0.5]
    else
        prim = [0.125, 0.0, 0.625]
    end

    u[i, ppp1, :] .= prim_conserve(prim, γ)
end

tspan = (0.0, 0.15)
prob = FREulerProblem(u, tspan, ps, γ, :dirichlet)
nt = tspan[2] ÷ dt |> Int
itg = init(prob, Midpoint(), saveat = tspan[2], adaptive = false, dt = dt)

@showprogress for iter = 1:nt
    @inbounds @threads for i = 1:ps.nx
        ũ = ps.iV * itg.u[i, :, 1]
        su = (ũ[end]^2) / (sum(ũ.^2) + 1e-6)
        isd = shock_detector(log10(su), ps.deg)
        λ = dt * exp(0.875/1 * (ps.deg)) * 0.15
        if isd
            for s = 1:3
                û = ps.iV * itg.u[i, :, s]
                FR.modal_filter!(û, λ; filter = :l2)
                #FR.modal_filter!(û, ℓ; filter = :lasso)
                #FR.modal_filter!(û, 4; filter = :exp)
                #FR.modal_filter!(û, 6; filter = :houli)
                itg.u[i, :, s] .= ps.V * û
            end
        end
    end

    step!(itg)
end

begin
    x = zeros(ncell * nsp)
    w = zeros(ncell * nsp, 3)
    for i = 1:ncell
        idx0 = (i - 1) * nsp

        for j = 1:nsp
            idx = idx0 + j
            x[idx] = ps.xpg[i, j]

            w[idx, :] .= itg.u[i, j, :]
        end
    end

    sol = zeros(ncell*nsp, 3)
    for i in axes(sol, 1)
        sol[i, :] .= conserve_prim(w[i, :], γ)
        sol[i, end] = 1 / sol[i, end]
    end

    plot(x, sol[:, 1], label="ρ", xlabel="x")
    plot!(x, sol[:, 2], label="u", xlabel="x")
    plot!(x, sol[:, 3], label="T", xlabel="x")
end


text:
function hS(x,Param)
f=0;
if Param.hS == "AgnesiAnnulus"
    f=Param.h/(((x-Param.phi0)/Param.aphi)^2+1);
elseif Param.hS == "AgnesiCart"
    f=Param.hC/(((x-Param.x0C)/Param.aC)^2+1);
else
    f=0;
end
return f
end


text:
"""
    rmkstruct(N, D; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (KRInfo, iz, nfp, ip)

Determine the Kronecker-structure and infinite pole-zero structure information of the rational matrix
`R(λ) := N(λ)./D(λ)` and return an `KRInfo` object, the multiplicities of infinite zeros in `iz`,
the number of finite poles in `nfp` and the multiplicities of infinite poles in `ip`.

The numerator `N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`,
for which the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

The denominator `D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`,
for which the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The information on the Kronecker-structure of the rational matrix `R(λ)` consists of
the right Kronecker indices `rki`, left Kronecker indices `lki`, the number of finite zeros `nf`
and the normal rank `nrank` and can be obtained from `KRInfo` as
`KRInfo.rki`, `KRInfo.lki`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`KRInfo`](@ref). Additionally, `KRInfo.id` contains the
infinite elementary divisors of the underlying system matrix pencil `S(λ)` used as linearization.

The Kronecker-structure information of the rational matrix `R(λ)` is obtained from the
Kronecker-structure information of the underlying linearization using the results of [1] and [2].
The determination of the Kronecker-structure information is performed by building an irreducible
descriptor system realization `(A-λE,B,C,D)` with `A-λE` a regular pencil of order `n`,
satisfying `R(λ) = C*inv(λE-A)*B+D` and reducing the structured system matrix pencil
`S(λ) = [A-λE B; C D]` to an appropriate Kronecker-like form (KLF) which exhibits
the number of its finite eigenvalues (also the finite zeros of `R(λ)`), the multiplicities of the
infinite eigenvalues (in excess with one to the multiplicities of infinite zeros of `R(λ)`),
the left and right Kronecker indices (also of `R(λ)`) and the normal rank
(in excess with `n` to the normal rank of `R(λ)`).
The Kronecker-structure information and pole-zero stucture information on `R(λ)` are recovered from
the Kronecker-structure information of `S(λ)` and `A-λE`, respectively, using the results of [1] and [2].

The right Kronecker indices are provided in the integer vector `rki`.
The number of elements of `rki` is the dimension of the right nullspace of the rational matrix `R(λ)`
and their sum is the least degree of a right polynomial nullspace basis or the least McMillan degree
of a right rational nullspace basis.

The left Kronecker indices are provided in the integer vector `lki`.
The number of elements of `lki` is the dimension of the left nullspace of the rational matrix `R(λ)`
and their sum is the least degree of a left polynomial nullspace basis or the least McMillan degree
of a left rational nullspace basis.

The multiplicities of infinite eigenvalues of `S(λ)` are provided in the integer vector `id`,
where each `i`-th element `id[i]` is the order of an infinite elementary divisor
(i.e., the multiplicity of an infinite eigenvalue). To each `id[i] > 1` corresponds an infinite zero
of `R(λ)` of multiplicity `id[i]-1`. The multiplicities of the infinite zeros of `R(λ)` are returned in `iz`.

The poles of `R(λ)` are the zeros of the regular pencil `A-λE`. The finite eigenvalues of `A-λE` are the
finite poles of `R(λ)` and their number is returned in `nfp`. The infinite zeros of `A-λE` are the
infinite poles of `R(λ)` and their  multiplicities are returned in `ip`.

The irreducible descriptor system based linearization is built using the methods described in [3] in conjunction with
pencil manipulation algorithms of [4] and [5] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The reductions to appropriate KLFs of `S(λ)` and `A-λE` are performed using orthogonal similarity transformations [5]
and involves rank decisions based  on rank revealing QR-decompositions with column pivoting, if `fast = true`,
or, the more reliable, SVD-decompositions, if `fast = false`.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances,  respectively, for the nonzero
coefficients of `N(λ)` and `D(λ)`.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `N(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `N(λ)`.

[1] G. Verghese, B. Levy, and T. Kailath, Generalized state-space systems, IEEE Trans. Automat. Control,
26:811-831 (1981).

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.

[3] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[4] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[5] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rmkstruct(N::AbstractArray{T1,3}, D::AbstractArray{T2,3}; fast::Bool = false, atol::Real = zero(real(T1)),
                   rtol::Real = (min(size(N)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
    A, E, B, C, D, blkdims = rm2ls(N, D; contr = true, obs = true, noseig = false,
                                   fast = fast, atol = atol, rtol = rtol)
    info = spkstruct(A, E, B, C, D; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
    iz  = info.id .- 1
    iz = iz[iz .> 0]
    n = size(A,1)
    info.nrank -= n
    nfp = blkdims[1]
    if n > nfp
       ii = nfp+1:n
       infop = pkstruct(view(A,ii,ii),view(E,ii,ii), fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
       ip = infop.id .- 1
       ip = ip[ip .> 0]
    else
      ip = Int[]
    end
    return info, iz, nfp, ip
end
function rmkstruct(P::AbstractArray{T,3}; kwargs...) where {T}
   info, iz, ip = pmkstruct(P; kwargs...)
   return info, iz, 0, ip
end
function rmkstruct(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
   D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return rmkstruct(poly2pm(N),poly2pm(D); kwargs...)
end
# function rmkstruct(N::Union{AbstractVecOrMat{<:Polynomial{T1,X}},Polynomial{T1,X},Number,AbstractVecOrMat{<:Number}},
#    D::Union{AbstractVecOrMat{<:Polynomial{T2,X}},Polynomial{T2,X},Number,AbstractVecOrMat{<:Number}}; kwargs...)  where {T1<:Number,T2<:Number,X}
#    return rmkstruct(poly2pm(N),poly2pm(D); kwargs...)
# end
# function rmkstruct(N::Union{AbstractVecOrMat{Polynomial{T1,X}}, AbstractVecOrMat{<:Polynomial}, Polynomial{T1,X},Polynomial,Number,AbstractVecOrMat{<:Number}},
#    D::Union{AbstractVecOrMat{Polynomial{T2,X}}, AbstractVecOrMat{<:Polynomial}, Polynomial{T2,X},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)  where {T1<:Number,T2<:Number,X}
#    return rmkstruct(poly2pm(N),poly2pm(D); kwargs...)
# end
function rmkstruct(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   info, iz, ip = pmkstruct(poly2pm(P); kwargs...)
   return info, iz, 0, ip
end
function rmkstruct(N::AbstractArray{T,3},
                   D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where T
   return rmkstruct(N,poly2pm(D); kwargs...)
end
function rmkstruct(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                   D::AbstractArray{T,3}; kwargs...) where T
   return rmkstruct(poly2pm(N),D; kwargs...)
end
"""
    rmzeros(N, D; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the finite and infinite zeros of the rational matrix `R(λ) := N(λ)./D(λ)` in `val`,
the multiplicities of infinite zeros in `iz` and the information on the Kronecker-structure of
the rational matrix `R(λ)` in the `KRInfo` object.

The information on the Kronecker-structure of the rational matrix `R(λ)` consists of
the right Kronecker indices `rki`, left Kronecker indices `lki`, the number of finite zeros `nf`
and the normal rank `nrank` and can be obtained from `KRInfo` as
`KRInfo.rki`, `KRInfo.lki`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`KRInfo`](@ref). Additionally, `KRInfo.id` contains the
infinite elementary divisors of the underlying linearization as a structured system matrix pencil.

The Kronecker-structure information of the rational matrix `R(λ)` is obtained from the
Kronecker-structure information of the underlying linearization using the results of [1] and [2].
The determination of the Kronecker-structure information is performed by building an irreducible
descriptor system realization `(A-λE,B,C,D)` with `A-λE` a regular pencil of order `n`,
satisfying `R(λ) = C*inv(λE-A)*B+D` and reducing the structured system matrix pencil
`S(λ) = [A-λE B; C D]` to an appropriate Kronecker-like form (KLF) which exhibits
the number of its finite eigenvalues (also the finite zeros of `R(λ)`), the multiplicities of the
infinite eigenvalues (in excess with one to the multiplicities of infinite zeros of `R(λ)`),
the left and right Kronecker indices (also of `R(λ)`) and the normal rank
(in excess with `n` to the normal rank of `R(λ)`).

The numerator `N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`,
for which the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

The denominator `D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`,
for which the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The irreducible descriptor system based linearization is built using the methods described in [3] in conjunction with
pencil manipulation algorithms of [4] and [5] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The reduction to the appropriate KLF of `S(λ)` is performed using orthogonal similarity transformations [5]
and involves rank decisions based on rank revealing QR-decompositions with column pivoting, if `fast = true`,
or, the more reliable, SVD-decompositions, if `fast = false`.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances,  respectively, for the nonzero
coefficients of `N(λ)` and `D(λ)`.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `N(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `N(λ)`.

[1] G. Verghese, B. Levy, and T. Kailath, Generalized state-space systems, IEEE Trans. Automat. Control,
26:811-831 (1981).

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.

[3] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[4] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[5] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rmzeros(N::AbstractArray{T1,3}, D::AbstractArray{T2,3};
                 fast = false, atol::Real = zero(real(T1)),
                 rtol::Real = (min(size(N)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
    A, E, B, C, D, blkdims = rm2ls(N, D; contr = true, obs = true, noseig = false,
                                   fast = fast, atol = atol, rtol = rtol)
    val, iz, info = spzeros(A, E, B, C, D; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
    info.nrank -= size(A,1)
    return val, iz, info
end
function rmzeros(N::AbstractArray{T,3}; kwargs...) where {T}
   return pmzeros2(N; kwargs...)
end
function rmzeros(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return rmzeros(poly2pm(N),poly2pm(D); kwargs...)
end
function rmzeros(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return pmzeros2(poly2pm(N); kwargs...)
end
function rmzeros(N::AbstractArray{T,3},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where T
   return rmzeros(N,poly2pm(D); kwargs...)
end
function rmzeros(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::AbstractArray{T,3}; kwargs...) where T
   return rmzeros(poly2pm(N),D; kwargs...)
end
"""
    rmpoles(N, D; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, ip, id)

Return the finite and infinite poles of the rational matrix `R(λ) := N(λ)./D(λ)`  in `val`,
the multiplicities of infinite poles in `ip` and the infinite elementary divisors of the pole pencil
of the underlying irreducible descriptor system based linearization of `R(λ)` in `id`.

The information on the finite and infinite poles of the rational matrix `R(λ)` is obtained from the
Kronecker-structure information of the underlying linearization using the results of [1] and [2].
The determination of the Kronecker-structure information is performed by building an irreducible
descriptor system realization `(A-λE,B,C,D)` with `A-λE` a regular pencil of order `n`,
satisfying `R(λ) = C*inv(λE-A)*B+D` and reducing the regular pencil `A-λE`
to an appropriate Kronecker-like form (KLF) which exhibits
the its finite eigenvalues (also the finite poles of `R(λ)`) and the multiplicities of the
infinite eigenvalues (in excess with one to the multiplicities of infinite poles of `R(λ)`).
The multiplicities of the infinite zeros of `A-λE` and of infinite poles of `R(λ)` are the same [2]
and are returned in `ip`. The number of infinite poles in `val` is equal to the sum of multiplicites in `ip`.
The multiplicities of the infinite eigenvalues of `A-λE` are returned in `id`.

The numerator `N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`,
for which the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

The denominator `D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`,
for which the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The irreducible descriptor system based linearization is built using the methods described in [3] in conjunction with
pencil manipulation algorithms of [4] and [5] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The reduction of `A-λE` to the KLF is performed using orthonal similarity transformations and involves rank decisions
based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances for the nonzero
coefficients of `N(λ)` and `D(λ)`,  respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `N(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `N(λ)`.

[1] G. Verghese, B. Levy, and T. Kailath, Generalized state-space systems, IEEE Trans. Automat. Control,
26:811-831 (1981).

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.

[3] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[4] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[5] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rmpoles(N::AbstractArray{T1,3}, D::AbstractArray{T2,3};
                 fast = false, atol::Real = zero(real(T1)),
                 rtol::Real = (min(size(N)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
   A, E, _, _, _, blkdims = rm2ls(N, D; contr = true, obs = true, noseig = false,
                                  fast = fast, atol = atol, rtol = rtol)
   n = size(A,1)
   nfp = blkdims[1]
   if nfp < n
      ii = nfp+1:n
      infop = pkstruct(view(A,ii,ii),view(E,ii,ii), fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      ip = infop.id .- 1
      ip = ip[ip .> 0]
      id = infop.id
   else
      ip = Int[]
      id = Int[]
   end
   return [eigvalsnosort!(view(A,1:nfp,1:nfp)); Inf*ones(eltype(A),sum(ip))], ip, id
end
function rmpoles(N::AbstractArray{T,3}; kwargs...) where T
   return pmpoles2(N; kwargs...)
end
function rmpoles(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return rmpoles(poly2pm(N),poly2pm(D); kwargs...)
end
function rmpoles(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return pmpoles2(poly2pm(N); kwargs...)
end
function rmpoles(N::AbstractArray{T,3},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where T
   return rmpoles(N,poly2pm(D); kwargs...)
end
function rmpoles(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::AbstractArray{T,3}; kwargs...) where T
   return rmpoles(poly2pm(N),D; kwargs...)
end
"""
    rmzeros1(N, D; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the finite and infinite zeros of the rational matrix `R(λ) := N(λ)./D(λ)` in `val`,
the multiplicities of infinite zeros in `iz` and the information on the Kronecker-structure of
the rational matrix `R(λ)` in the `KRInfo` object.

The information on the Kronecker-structure of the rational matrix `R(λ)` consists of
the right Kronecker indices `rki`, left Kronecker indices `lki`, the number of finite zeros `nf`
and the normal rank `nrank` and can be obtained from `KRInfo` as
`KRInfo.rki`, `KRInfo.lki`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`KRInfo`](@ref). Additionally, `KRInfo.id` contains the
infinite elementary divisors of the underlying pencil based linearization.

The Kronecker-structure information of the rational matrix `R(λ)` is obtained from the
Kronecker-structure information of the underlying pencil based linearization using the results of [1] and [2].
The determination of the Kronecker-structure information is performed by building a strongly minimal
pencil based linearization `(A-λE,B-λF,C-λG,D-λH)` with `A-λE` a regular pencil of order `n`,
satisfying `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH` and reducing the structured system matrix pencil
`S(λ) = [A-λE B-λF; C-λG D-λH]` to an appropriate Kronecker-like form (KLF) which exhibits
the number of its finite eigenvalues (also the finite zeros of `R(λ)`), the multiplicities of the
infinite eigenvalues (in excess with one to the multiplicities of infinite zeros of `R(λ)`),
the left and right Kronecker indices (also of `R(λ)`) and the normal rank
(in excess with `n` to the normal rank of `R(λ)`).

The numerator `N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`,
for which the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

The denominator `D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`,
for which the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The strongly minimal pencil based linearization is built using the methods described in [3] in conjunction with
pencil manipulation algorithms of [4] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The reduction to the appropriate KLF of `S(λ)` is performed using orthogonal similarity transformations [5]
and involves rank decisions based on rank revealing QR-decompositions with column pivoting, if `fast = true`,
or, the more reliable, SVD-decompositions, if `fast = false`.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances,  respectively, for the nonzero
coefficients of `N(λ)` and `D(λ)`.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `N(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `N(λ)`.

[1] G. Verghese, B. Levy, and T. Kailath, Generalized state-space systems, IEEE Trans. Automat. Control,
26:811-831 (1981).

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.

[3] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[4] F.M. Dopico, M.C. Quintana and P. Van Dooren, Linear system matrices of rational transfer functions, to appear in "Realization and Model Reduction of Dynamical Systems",
A Festschrift to honor the 70th birthday of Thanos Antoulas", Springer-Verlag. [arXiv:1903.05016](https://arxiv.org/pdf/1903.05016.pdf)

[5] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rmzeros1(N::AbstractArray{T1,3}, D::AbstractArray{T2,3};
                  fast = false, atol::Real = zero(real(T1)),
                   rtol::Real = (min(size(N)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
   A, E, B, F, C, G, D, H, blkdims = rm2lps(N, D; obs = true, contr = true, fast = fast,
                                            atol = atol, rtol = rtol)
   val, iz, info = pzeros([A B; C D], [E F; G H]; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   info.nrank -= size(A,1)
   return val, iz, info
end
function rmzeros1(N::AbstractArray{T,3}; kwargs...) where {T}
   return pmzeros1(N; kwargs...)
end
function rmzeros1(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return rmzeros1(poly2pm(N),poly2pm(D); kwargs...)
end
function rmzeros1(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return pmzeros1(poly2pm(N); kwargs...)
end
function rmzeros1(N::AbstractArray{T,3},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where T
   return rmzeros1(N,poly2pm(D); kwargs...)
end
function rmzeros1(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::AbstractArray{T,3}; kwargs...) where T
   return rmzeros1(poly2pm(N),D; kwargs...)
end
"""
    rmpoles1(N, D;  fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, ip, id)

Return the finite and infinite poles of the rational matrix `R(λ) := N(λ)./D(λ)`  in `val`,
the multiplicities of infinite poles in `ip` and the infinite elementary divisors of the pole pencil
of the underlying strongly irreducible pencil based linearization of  `R(λ)` in `id`.

The information on the finite and infinite poles of the rational matrix `R(λ)` is obtained from the
pole-structure information of the underlying pencil based linearization using the results of [1] and [2].
The determination of the pole-structure information is performed by building a strongly minimal
pencil based linearization `(A-λE,B-λF,C-λG,D-λH)` with `A-λE` a regular pencil of order `n`,
satisfying `R(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH` and reducing the structured system matrix pencil
`Sp(λ) = [A-λE -λF; -λG -λH]` to an appropriate Kronecker-like form (KLF) which exhibits
the number of its finite eigenvalues (also the finite poles of `R(λ)`) and the multiplicities of the
infinite eigenvalues (in excess with one to the multiplicities of infinite poles of `R(λ)`).
The multiplicities of the infinite zeros of `Sp(λ)` and of infinite poles of `R(λ)` are the same [2]
and are returned in `ip`. The number of infinite poles in `val` is equal to the sum of multiplicites in `ip`.
The multiplicities of the infinite eigenvalues of `Sp(λ)` are returned in `id`.

The numerator `N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`,
for which the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

The denominator `D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`,
for which the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The strongly minimal pencil based linearization is built using the methods described in [3] in conjunction with
pencil manipulation algorithms of [4] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The reduction of `Sp(λ)` to the KLF is performed using orthonal similarity transformations [5] and involves rank decisions
based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances,  respectively, for the nonzero
coefficients of `N(λ)` and `D(λ)`.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `N(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `N(λ)`.

[1] G. Verghese, B. Levy, and T. Kailath, Generalized state-space systems, IEEE Trans. Automat. Control,
26:811-831 (1981).

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.

[3] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[4] F.M. Dopico, M.C. Quintana and P. Van Dooren, Linear system matrices of rational transfer functions, to appear in "Realization and Model Reduction of Dynamical Systems",
A Festschrift to honor the 70th birthday of Thanos Antoulas", Springer-Verlag. [arXiv:1903.05016](https://arxiv.org/pdf/1903.05016.pdf)

[5] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rmpoles1(N::AbstractArray{T1,3}, D::AbstractArray{T2,3};
                  fast = false, atol::Real = zero(real(T1)),
                  rtol::Real = (min(size(N)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
   A, E, _, F, _, G, _, H, blkdims = rm2lps(N,D; obs = true, contr = true, fast = fast, atol = atol, rtol = rtol)
   n = size(A,1)
   nfp = blkdims[1]
   p, m = size(H)
   ii = nfp+1:n
   ni = n-nfp
   T = promote_type(T1,T2)
   M = [A[ii,ii] zeros(T,ni,p+m); zeros(T,p,ni+m) I; zeros(T,m,ni) I zeros(T,m,p)]
   N = [E[ii,ii] F[ii,:] zeros(T,ni,p); G[:,ii] H zeros(T,p,p); zeros(T,m,ni+m+p)]
   infop = pkstruct(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   ip = infop.id .- 1
   ip = ip[ip .> 0]
   id = infop.id
   return [eigvalsnosort!(view(A,1:nfp,1:nfp)); Inf*ones(eltype(A),sum(ip))], ip, id
   # p, m = size(D)
   # n = size(A,1)
   # T = eltype(A)
   # M = [A zeros(T,n,p+m); zeros(T,p,n+m) I; zeros(T,m,n) I zeros(T,m,p)]
   # N = [E F zeros(T,n,p); G H zeros(T,p,p); zeros(T,m,n+m+p)]
   # val, iz, info = pzeros(M, N; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   # return val, iz, info.id
end
function rmpoles1(N::AbstractArray{T,3}; kwargs...) where T
   return pmpoles1(N; kwargs...)
end
function rmpoles1(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where {T1,T2}
   return rmpoles1(poly2pm(N),poly2pm(D); kwargs...)
end
function rmpoles1(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where {T}
   return pmpoles1(poly2pm(N); kwargs...)
end
function rmpoles1(N::AbstractArray{T,3},
                 D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where T
   return rmpoles1(N,poly2pm(D); kwargs...)
end
function rmpoles1(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                 D::AbstractArray{T,3}; kwargs...) where T
   return rmpoles1(poly2pm(N),D; kwargs...)
end
"""
       rmrank(N, D; fastrank = true, atol = 0, rtol = atol > 0 ? 0 : n*ϵ)

Determine the normal rank of the rational matrix `R(λ) := N(λ)./D(λ)`.

If `fastrank = true`, the rank is evaluated by counting how many singular values of `R(γ)` have magnitude greater
than `max(atol, rtol*σ₁)`, where `σ₁` is the largest singular value of `R(γ)` and `γ` is a randomly generated
complex value of magnitude equal to one.
If `fastrank = false`, first a structured linearization of `R(λ)` is built in the form of a descriptor system matrix
`S(λ) = [A-λE B; C D]` with `A-λE` an order `n` regular subpencil, and then its normal rank `nr` is determined.
The normal rank of `R(λ)` is `nr - n`.

The numerator `N(λ)` is a polynomial matrix of the form `N(λ) = N_1 + λ N_2 + ... + λ**k N_(k+1)`,
for which the coefficient matrices `N_i`, `i = 1, ..., k+1` are stored in the 3-dimensional matrix `N`,
where `N[:,:,i]` contains the `i`-th coefficient matrix `N_i` (multiplying `λ**(i-1)`).

The denominator `D(λ)` is a polynomial matrix of the form `D(λ) = D_1 + λ D_2 + ... + λ**l D_(l+1)`,
for which the coefficient matrices `D_i`, `i = 1, ..., l+1`, are stored in the 3-dimensional matrix `D`,
where `D[:,:,i]` contain the `i`-th coefficient matrix `D_i` (multiplying `λ**(i-1)`).

Alternatively, `N(λ)` and `D(λ)` can be specified as matrices of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
In this case, no check is performed that `N(λ)` and `D(λ)` have the same indeterminates.

The irreducible descriptor system based linearization is built using the methods described in [1] in conjunction with
pencil manipulation algorithms of [2] and [3] to compute reduced order linearizations. These algorithms
employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition, if `fast = false`.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

The reduction to the appropriate KLF of `S(λ)` is performed using orthogonal similarity transformations [3]
and involves rank decisions based on rank revealing SVD-decompositions.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances for the nonzero
coefficients of `N(λ)` and `D(λ)`,  respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `N(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `N(λ)`.

[1] A. Varga, On computing the Kronecker structure of polynomial and rational matrices using Julia, 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).

[2] P. Van Dooreen, The generalized eigenstructure problem in linear system theory,
IEEE Transactions on Automatic Control, vol. AC-26, pp. 111-129, 1981.

[3] A. Varga, Solving Fault Diagnosis Problems - Linear Synthesis Techniques, Springer Verlag, 2017.
"""
function rmrank(N::AbstractArray{T1,3}, D::AbstractArray{T2,3}; fastrank::Bool = true, atol::Real = zero(real(T1)),
                rtol::Real = (min(size(N)...)*eps(real(float(one(T1)))))*iszero(atol)) where {T1,T2}
   if fastrank
       return rank(rmeval(N,D,exp(rand()*im)), atol = atol, rtol = rtol)
   else
      A, E, B, C, D, blkdims = rm2ls(N, D; fast = false, atol = atol, rtol = rtol)
      r = sprank(A, E, B, C, D, fastrank = false, atol1 = atol, atol2 = atol, rtol = rtol)
      return r-size(A,1)
   end
end
function rmrank(N::AbstractArray{T,3}; kwargs...) where T
    pmrank(N; kwargs...)
end
function rmrank(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
   return rmrank(poly2pm(N),poly2pm(D); kwargs...)
end
function rmrank(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...)
    pmrank(poly2pm(N); kwargs...)
end
function rmrank(N::AbstractArray{T,3},
                D::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where T
   return rmrank(N,poly2pm(D); kwargs...)
end
function rmrank(N::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}},
                D::AbstractArray{T,3}; kwargs...) where T
   return rmrank(poly2pm(N),D; kwargs...)
end