text:
const DimOrDimTuple = Union{Dimension,Tuple{Vararg{<:Dimension}}}
const IntOrIntTuple = Union{Int,Tuple{Vararg{<:Int}}}

"""
    aggregate(method, object, scale)

Aggregate array, or all arrays in a stack or series, by some scale.

# Arguments

- `method`: a function such as `mean` or `sum` that can combine the
  value of multiple cells to generate the aggregated cell, or a [`Locus`]($DDlocusdocs)
  like `Start()` or `Center()` that species where to sample from in the interval.
- `object`: Object to aggregate, like `AbstractGeoSeries`, `AbstractStack`,
  `AbstractGeoArray`, `Dimension`
- `scale`: the aggregation factor, which can be an integer, a tuple of integers
  for each dimension, or any `Dimension`, `Selector` or `Int` combination you can
  usually use in `getindex`. Using a `Selector` will determine the scale by the
  distance from the start of the index.
"""
function aggregate end

"""
    aggregate(method, series::AbstractGeoSeries, scale)

Aggregate an [`AbstractGeoSeries`](@ref) by `scale` using `method`.

Returns a [`GeoSeries`](ref).
"""
function aggregate(
    method, series::AbstractGeoSeries, scale, args...; progress=true, kwargs...
)
    f = i -> aggregate(method, series[i], scale, args...; progress=false, kwargs...)
    data = if progress
        @showprogress "Aggregating series..." map(f, 1:length(series))
    else
        map(f, 1:length(series))
    end
    return rebuild(series, data)
end
"""
    aggregate(method, stack::AbstractGeoStack, scale)

Aggregate an [`AbstractGeoStack`](@ref) by `scale` using `method`.

Returns a [`GeoStack`](ref).
"""
function aggregate(
    method, stack::AbstractGeoStack, scale; keys=keys(stack), progress=true
)
    f = key -> aggregate(method, stack[key], scale)
    keys_nt = NamedTuple{keys}(keys)
    data = if progress
        @showprogress "Aggregating stack..." map(f, keys_nt)
    else
        map(f, keys_nt)
    end
    return GeoStack(stack; data=data)
end

# DimensionalData methods

"""
    aggregate(method, src::AbstractDimArray, scale)

Aggregate an `AbstractDimArray` by `scale` using `method`.

[`DiskGeoArray`](@ref) will be converted to [`GeoArray`](@ref).
"""
aggregate(method, src::AbstractDimArray, scale) =
    aggregate!(method, alloc_ag(method, src, scale), src, scale)
aggregate(method, src::DiskGeoArray, scale) =
    aggregate(method, GeoArray(src), scale)

"""
    aggregate(method, dim::Dimension, scale)

Aggregate a `Dimension` by `scale` using `method`.
"""
function aggregate(method, dim::Dimension, scale)
    intscale = scale2int(dim, scale)
    intscale == 1 && return dim
    start, stop = endpoints(dim, method, intscale)
    return rebuild(dim, val(dim)[start:scale:stop], aggregate(method, mode(dim), intscale))
end

function endpoints(dim, method, scale)
    start = firstindex(dim) + agoffset(method, dim, scale)
    stop = (length(dim) ÷ scale) * scale
    return start, stop
end

"""
    aggregate(method, dim::IndexMode, scale)

Aggregate an `IndexMode` by `scale` using `method`.
"""
aggregate(method, mode::IndexMode, scale) = mode
aggregate(method, mode::AbstractSampled, scale) =
    rebuild(mode; span=aggregate(method, span(mode), scale))

"""
    aggregate(method, dim::Span, scale)

Aggregate a `Span` by `scale` using `method`.
"""
aggregate(method, span::Span, scale) = span
aggregate(method, span::Regular, scale) = Regular(val(span) * scale)

"""
    aggregate!(method, dst::AbstractDimArray, src::AbstractDimArray, scale)

Aggregate array `src` to array `dst` by `scale`, using `method`.

# Arguments

- `method`: a function such as `mean` or `sum` that can combine the
  value of multiple cells to generate the aggregated cell, or a [`Locus`]($DDlocusdocs)
  like `Start()` or `Center()` that species where to sample from in the interval.
- `scale`: the aggregation factor, which can be an integer, a tuple of integers
  for each dimension, or any `Dimension`, `Selector` or `Int` combination you can
  usually use in `getindex`. Using a `Selector` will determine the scale by the
  distance from the start of the index in the `src` array.
"""
aggregate!(locus::Locus, dst::AbstractDimArray, src, scale) =
    aggregate!((locus,), dst, src, scale)
function aggregate!(loci::Tuple{Locus,Vararg}, dst::AbstractDimArray, src, scale)
    intscale = scale2int(dims(src), scale)
    offsets = agoffset.(loci, intscale)
    for I in CartesianIndices(dst)
        dst[I] = src[(upsample.(Tuple(I), intscale) .+ offsets)...]
    end
    return dst
end
# Function/functor methods
function aggregate!(f, dst::AbstractDimArray, src, scale)
    intscale = scale2int(dims(src), scale)
    for I in CartesianIndices(dst)
        topleft = upsample.(Tuple(I), intscale)
        bottomright = topleft .+ intscale .- 1
        dst[I] = f(view(src, map(:, topleft, bottomright)...))
    end
    return dst
end

"""
    alloc_ag(method, A::AbstractDimArray, scale)

Allocate an array of the correct size to aggregate `A` by `scale`
"""
alloc_ag(method, A::AbstractDimArray, scale) = alloc_ag((method,), A, scale)
function alloc_ag(method::Tuple, A::AbstractDimArray, scale)
    intscale = scale2int(dims(A), scale)
    # Aggregate the dimensions
    dims_ = aggregate.(method, dims(A), intscale)
    # Dim aggregation determines the array size
    data_ = similar(data(A), map(length, dims_)...)
    return rebuild(A; data=data_, dims=dims_)
end

# Convert scale or tuple of scale to integer using dims2indices
scale2int(dims, scale::Tuple) = dims2indices(dims, scale)
scale2int(dims, scale::Int) = scale

"""
    disaggregate(method, object, scale)

Disaggregate array, or all arrays in a stack or series, by some scale.

# Arguments

- `method`: a function such as `mean` or `sum` that can combine the
  value of multiple cells to generate the aggregated cell, or a [`Locus`]($DDlocusdocs)
  like `Start()` or `Center()` that species where to sample from in the interval.
- `object`: Object to aggregate, like `AbstractGeoSeries`, `AbstractStack`,
  `AbstractGeoArray`, `Dimension`
- `scale`: the aggregation factor, which can be an integer, a tuple of integers
  for each dimension, or any `Dimension`, `Selector` or `Int` combination you can
  usually use in `getindex`. Using a `Selector` will determine the scale by the
  distance from the start of the index.
"""
function disaggregate end

"""
    disaggregate(method, series::AbstractGeoSeries, scale)

Disaggregate an [`AbstractGeoSeries`](@ref) by `scale` using `method`.
"""
function disaggregate(method, series::AbstractGeoSeries, scale; progress=true, kwargs...)
    f = i -> disaggregate(method, series[i], scale; progress=false, kwargs...)
    return if progress
        @showprogress "Disaggregating series..." map(f, 1:length(series))
    else
        map(f, 1:length(series))
    end
end
"""
    disaggregate(method, stack::AbstractGeoStack, scale)

Disaggregate an [`AbstractGeoStack`](@ref) by `scale` using `method`.
"""
function disaggregate(method, stack::AbstractGeoStack, scale; keys=keys(stack), progress=true)
    f = key -> disaggregate(method, stack[key], scale)
    keys_nt = NamedTuple{keys}(keys)
    data = if progress
        @showprogress "Disaggregating stack..." map(f, keys_nt)
    else
        map(f, keys_nt)
    end
    return GeoStack(stack; data=data)
end

# DimensionalData methods
"""
    disaggregate(method, src::AbstractDimArray, scale)

Disaggregate an `AbstractDimArray` by `scale` using `method`.

[`DiskGeoArray`](@ref) will be converted to [`GeoArray`](@ref).
"""
disaggregate(method, src::AbstractDimArray, scale) =
    disaggregate!(method, alloc_disag(method, src, scale), src, scale)
disaggregate(method, src::DiskGeoArray, scale) = disaggregate(method, GeoArray(src), scale)

"""
    disaggregate(method, dim::Dimension, scale)

Disaggregate a `Dimension` by `scale` using `method`.
"""
function disaggregate(locus::Locus, dim::Dimension, scale)
    intscale = scale2int(dim, scale)
    intscale == 1 && return dim
    len = length(dim) * intscale
    step_ = step(mode(dim)) / intscale
    start = dim[1] - agoffset(locus, intscale) * step_
    stop = start + (len - 1)  * step_
    index = LinRange(start, stop, len)
    return rebuild(dim, index, disaggregate(locus, mode(dim), intscale))
end
disag_index(locus::Start, dim, scale) =
    LinRange(dim[1], dim[end] + (scale - 1) * step(dim), length(dim) * scale)
disag_index(locus::End, dim, scale) =
    LinRange(dim[1], dim[end] + (scale - agoffset(locus)) * step(dim), length(dim) * scale)

"""
    disaggregate(method, dim::IndexMode, scale)

Disaggregate an `IndexMode` by `scale` using `method`.
"""
disaggregate(method, mode::IndexMode, scale) = mode
disaggregate(method, mode::AbstractSampled, scale) =
    rebuild(mode; span=disaggregate(method, span(mode), scale))

"""
    disaggregate(method, dim::Span, scale)

Disaggregate a `Span` by `scale` using `method`.
"""
disaggregate(method, span::Span, scale) = span
disaggregate(method, span::Regular, scale) = Regular(val(span) / scale)

"""
    disaggregate!(method, dst::AbstractDimArray, src::AbstractDimArray, scale)

Disaggregate array `src` to array `dst` by some scale, using `method`.

- `method`: a function such as `mean` or `sum` that can combine the
  value of multiple cells to generate the aggregated cell, or a [`Locus`]($DDlocusdocs)
  like `Start()` or `Center()` that species where to sample from in the interval.
- `scale`: the aggregation factor, which can be an integer, a tuple of integers
  for each dimension, or any `Dimension`, `Selector` or `Int` combination you can
  usually use in `getindex`. Using a `Selector` will determine the scale by the
  distance from the start of the index in the `src` array.
"""
disaggregate!(locus::Locus, dst::AbstractDimArray, src, scale) =
    disaggregate!((locus,), dst, src, scale)
function disaggregate!(loci::Tuple{Locus,Vararg}, dst::AbstractDimArray, src, scale)
    intscale = scale2int(dims(src), scale)
    for I in CartesianIndices(dst)
        dst[I] = src[(downsample.(Tuple(I), intscale))...]
    end
    return dst
end

"""
    alloc_ag(method, A::AbstractDimArray, scale)

Allocate an array of the correct size to disaggregate `A` by `scale`
"""
alloc_disag(method, A::AbstractDimArray, scale) = alloc_disag((method,), A, scale)
function alloc_disag(method::Tuple, A::AbstractDimArray, scale)
    intscale = scale2int(dims(A), scale)
    dims_ = disaggregate.(method, dims(A), intscale)
    # Dim aggregation determines the array size
    data_ = similar(data(A), map(length, dims_)...)
    return rebuild(A; data=data_, dims=dims_)
end

"""
    upsample(index::Int, scale::Int)

Convert indicies from the aggregated array to the larger original array.
"""
upsample(index::Int, scale::Int) = (index - 1) * scale + 1

"""
    downsample(index::Int, scale::Int)

Convert indicies from the original array to the aggregated array.
"""
downsample(index::Int, scale::Int) = (index - 1) ÷ scale + 1

agoffset(locus::Locus, dim::Dimension, scale) = agoffset(locus, scale)
agoffset(method, dim::Dimension, scale) = agoffset(locus(dim), scale)
agoffset(locus::Start, scale) = 0
agoffset(locus::End, scale) = scale - 1
agoffset(locus::Center, scale) = scale ÷ 2


text:
module HypergeoMat

export hypergeomPQ
export lmvgamma
export mvgamma
export BesselA

module HypergeomPQ

export hypergeomPQ
import LinearAlgebra

function DictParts(m::Integer, n::Integer)
  D = Dict{Int64,Int64}()
  Last = hcat([0, m, m])
  fin = 0
  for i = 1:n
    NewLast = Array{Int64}(undef, 3, 0)
    for record in eachcol(Last)
      @inbounds manque = record[2]
      @inbounds l = min(manque, record[3])
      if l > 0
        @inbounds D[record[1]+1] = fin + 1
        x = Array{Int64}(undef, 3, l)
        for j = 1:l
          @inbounds x[:, j] = [fin + j, manque - j, j]
        end
        NewLast = hcat(NewLast, x)
        fin += l
      end
    end
    Last = NewLast
  end
  return (D, fin)
end

function _N(dico::Dict{I,I}, lambda::Vector{I}) where {I<:Integer}
  lambda = filter(x -> x > 0, lambda)
  if isempty(lambda)
    return 1
  end
  @inbounds dico[_N(dico, lambda[1:end-1])] + last(lambda)
end

function _T(
  alpha::Real,
  a::Vector{<:Number},
  b::Vector{<:Number},
  kappa::Vector{<:Integer}
)
  if isempty(kappa) || kappa[1] == 0
    return 1
  end
  i = lastindex(kappa)
  @inbounds c = kappa[i] - 1 - (i - 1) / alpha
  @inbounds d = alpha * kappa[i] - i
  @inbounds s = collect(1:(kappa[i]-1))
  e = d .- alpha .* s + map(j -> count(>=(j), kappa), s)
  g = e .+ 1
  ss = 1:(i-1)
  @inbounds f = alpha .* kappa[ss] - collect(ss) .- d
  h = f .+ alpha
  l = h .* f
  prod1 = prod(a .+ c)
  prod2 = prod(b .+ c)
  prod3 = prod((g .- alpha) .* e ./ (g .* (e .+ alpha)))
  prod4 = prod((l - f) ./ (l + h))
  out = prod1 / prod2 * prod3 * prod4
  return isinf(out) || isnan(out) ? 0 : out
end

function betaratio(
  kappa::Vector{I},
  mu::Vector{I},
  k::I,
  alpha::T,
) where {I<:Integer,T<:Real}
  @inbounds muk = mu[k]
  t = k - alpha * muk
  @inbounds u = map(i -> t + 1 - i + alpha * kappa[i], 1:k)
  @inbounds v = map(i -> t - i + alpha * mu[i], 1:(k-1))
  muPrime = dualPartition(mu)
  @inbounds w = map(i -> muPrime[i] - t - alpha * i, 1:(muk-1))
  prod1 = prod(u ./ (u .+ alpha .- 1))
  prod2 = prod((v .+ alpha) ./ v)
  prod3 = prod((w .+ alpha) ./ w)
  return alpha * prod1 * prod2 * prod3
end

function dualPartition(lambda::Vector{T}) where {T<:Integer}
  out = T[]
  if !isempty(lambda)
    @inbounds for i in 1:lambda[1]
      push!(out, count(>=(i), lambda))
    end
  end
  return out
end

# ------------------------------------------------------------------------------
function hypergeomI(
  m::I,
  alpha::Real,
  a::Vector{<:Number},
  b::Vector{<:Number},
  n::I,
  x::Number,
)::Number where {I<:Integer}
  function summation(
    i::Int64,
    z::Number,
    j::I,
    kappa::Vector{I},
  )
    function go(kappai::I, zz::Number, s::Number)
      @inbounds if i == 0 && kappai > j || i > 0 && kappai > min(kappa[i], j)
        return s
      else
        kappap = vcat(kappa, [kappai])
        t = _T(alpha, a, b, filter(x -> x > 0, kappap))
        zp = zz * x * (n - i + alpha * (kappai - 1)) * t
        sp = j > kappai && i <= n ? s + summation(i + 1, zp, j - kappai, kappap) : s
        spp = sp + zp
        go(kappai + 1, zp, spp)
      end
    end # end go ---------------------------------------------------------------
    go(I(1), z, 0)
  end # end summation ----------------------------------------------------------
  return 1 + summation(0, 1, m, I[])
end

# ------------------------------------------------------------------------------
"""
    hypergeomPQ(m, a, b, x[, alpha])

Compute the truncated hypergeometric function of a matrix argument given the
  eigen values of the matrix.

# Arguments
- `m`: truncation weight of the summation, a positive integer
- `a`: the "upper" parameters, a real or complex vector, possibly empty
- `b`: the "lower" parameters, a real or complex vector, possibly empty
- `x`: real or complex vector, the eigen values
- `alpha`: the alpha parameter, a positive number; if missing, `alpha=2` is used
"""
function hypergeomPQ(
  m::Integer,
  a::Vector{<:Number},
  b::Vector{<:Number},
  x::Vector{<:Number},
  alpha::Union{Nothing,Real} = nothing,
)
  if isnothing(alpha)
    alpha = 2
  end
  n = length(x)
  @inbounds if all(x .== x[1])
    @inbounds return hypergeomI(m, alpha, a, b, n, x[1])
  end
  function jack(
    k::I,
    beta::Number,
    c::I,
    t::I,
    mu::Vector{I},
    jarray::Array{<:Number,2},
    kappa::Vector{I},
    nkappa::I,
  ) where {I<:Integer}
    lmu = length(mu)
    for i in ((max(k, 1):count(>(0), mu)))
      @inbounds u = mu[i]
      @inbounds if lmu == i || u > mu[i+1]
        gamma = beta * betaratio(kappa, mu, i, alpha)
        mup = copy(mu)
        @inbounds mup[i] = u - 1
        filter!(>(0), mup)
        if length(mup) >= i && u > 1
          jarray = jack(i, gamma, c + 1, t, mup, jarray, kappa, nkappa)
        else
          if nkappa > 1
            if !isempty(mup)
              @inbounds jarray[nkappa, t] +=
                gamma * jarray[_N(dico, mup)-1, t-1] * x[t]^(c + 1)
            else
              @inbounds jarray[nkappa, t] += gamma * x[t]^(c + 1)
            end
          end
        end
      end
    end
    if k == 0
      if nkappa > 1
        @inbounds jarray[nkappa, t] += jarray[nkappa, t-1]
      end
    else
      @inbounds jarray[nkappa, t] += beta * x[t]^c * jarray[_N(dico, mu)-1, t-1]
    end
    return jarray
  end # end jack ---------------------------------------------------------------
  function summation(
    i::I,
    z::Number,
    j::I,
    kappa::Vector{I},
    jarray::Array{<:Number,2},
  ) where {I<:Integer}
    function go(kappai::I, zp::Number, s::Number)
      if i == n || i == 0 && kappai > j || i > 0 && kappai > min(last(kappa), j)
        return s
      end
      kappap = vcat(kappa, [kappai])
      nkappa = _N(dico, kappap) - 1
      zpp = zp * _T(alpha, a, b, kappap)
      @inbounds if nkappa > 1 && (length(kappap) == 1 || kappap[2] == 0)
        @inbounds jarray[nkappa, 1] =
          x[1] * (1 + alpha * (kappap[1] - 1)) * jarray[nkappa-1, 1]
      end
      for t in 2:n
        jarray = jack(I(0), 1, I(0), t, kappap, jarray, kappap, nkappa)
      end
      @inbounds sp = s + zpp * jarray[nkappa, n]
      if j > kappai && i <= n
        spp = summation(i + 1, zpp, j - kappai, kappap, jarray)
        go(kappai + 1, zpp, sp + spp)
      else
        go(kappai + 1, zpp, sp)
      end
    end # end go ---------------------------------------------------------------
    go(I(1), z, 0)
  end # end summation ----------------------------------------------------------
  (dico, Pmn) = DictParts(m, n)
  T = typeof(x[1])
  J = zeros(T, Pmn, n)
  @inbounds J[1, :] = cumsum(x)
  return 1 + summation(0, T(1), m, Int64[], J)
end

# univariate -------------------------------------------------------------------
"""
    hypergeomPQ(m, a, b, x)

Compute the truncated hypergeometric function of a scalar argument.

# Arguments
- `m`: truncation weight of the summation, a positive integer
- `a`: the "upper" parameters, a real or complex vector, possibly empty
- `b`: the "lower" parameters, a real or complex vector, possibly empty
- `x`: scalar, real or complex
"""
function hypergeomPQ(
  m::Integer,
  a::Vector{<:Number},
  b::Vector{<:Number},
  x::Number
)
  return hypergeomI(m, 2, a, b, 1, x)
end

# matrix argument --------------------------------------------------------------
"""
    hypergeomPQ(m, a, b, X[, alpha])

Compute the truncated hypergeometric function of a matrix argument. The
  hypergeometric function is usually defined for a symmetric real matrix
  or a Hermitian complex matrix but arbitrary square matrices are allowed.

# Arguments
- `m`: truncation weight of the summation, a positive integer
- `a`: the "upper" parameters, a real or complex vector, possibly empty
- `b`: the "lower" parameters, a real or complex vector, possibly empty
- `X`: a square matrix, real or complex
- `alpha`: the alpha parameter, a positive number; if missing, `alpha=2` is used
"""
function hypergeomPQ(
  m::Integer,
  a::Vector{<:Number},
  b::Vector{<:Number},
  X::Matrix{<:Number},
  alpha::Union{Nothing,Real} = nothing
)
  x = LinearAlgebra.eigvals(X)
  return hypergeomPQ(m, a, b, x, alpha)
end

end # end module HypergeomPQ

module Mvgamma

export lmvgamma
export mvgamma
import GSL

"""
    lmvgamma(z, p)

Compute the logarithm of the multivariate Gamma function.

# Arguments
- `z`: real or complex number
- `p`: positive integer, the dimension
"""
function lmvgamma(
  z::Union{R,T},
  p::Integer
) where {R<:Real,T<:Complex{R}}
  if real(z) <= 0
    throw(DomainError(z, "The real part of `z` is not positive."))
  end
  C = p*(p-1)/4.0 * log(pi)
  isComplex = eltype(z) <: Complex
  if isComplex
    z_re = real(z)
    z_im = imag(z)
    (a, b) = GSL.sf_lngamma_complex_e(z_re, z_im)
    S = complex(a.val, b.val)
    for i in 2:p
      (a, b) = GSL.sf_lngamma_complex_e(z_re + (1 - i)/2.0, z_im)
      S = S + complex(a.val, b.val)
    end
  else
    S = GSL.sf_lngamma(z)
    for i in 2:p
      S = S + GSL.sf_lngamma(z + (1 - i)/2.0)
    end
  end
  return C + S
end

"""
    mvgamma(z, p)

Compute the multivariate Gamma function.

# Arguments
- `z`: real or complex number
- `p`: positive integer, the dimension
"""
function mvgamma(
  z::Union{R,T},
  p::Integer
) where {R<:Real,T<:Complex{R}}
  if imag(z) == 0 && real(z) <= 0 && trunc(z) == z
    throw(DomainError(z, "The argument `z` is a non-positive integer."))
  end
  if real(z) > 0.0
    result = exp(lmvgamma(z, p))
  else
    n = 1 + floor(-real(z))
    pochhammer = z
    for i in 2:n
      pochhammer = pochhammer * (z + i - 1.0)
    end
    result = exp(lmvgamma(z + n, p)) / pochhammer
  end
  return result
end

end # end module Mvgamma

module Bessel

export BesselA

"""
    BesselA(m, X, nu)

Compute the truncated Herz's type one Bessel function of a matrix argument. It
  is usually defined for a symmetric real matrix or a Hermitian complex matrix
  but arbitrary square matrices are allowed.

# Arguments
- `m`: truncation weight of the hypergeometric function, a positive integer
- `X`: a square matrix, real or complex
- `nu`: the order parameter, real or complex number with `real(nu)>-1`
"""
function BesselA(
  m::Integer,
  X::Matrix{Number},
  nu::Number
)
  if real(nu) <= -1
    throw(DomainError(nu, "The real part of `nu` is smaller than -1."))
  end
  p = size(X)[1]
  b = nu + (p+1) / 2
  return hypergeomPQ(m, Int64[], [b], -X, 2) / mvgamma(b, p)
end

"""
    BesselA(m, x, nu)

Compute the truncated Herz's type one Bessel function of a matrix argument
  given the eigen values of the matrix.

# Arguments
- `m`: truncation weight of the hypergeometric function, a positive integer
- `x`: the eigen values, a vector of real or complex numbers
- `nu`: the order parameter, real or complex number with `real(nu)>-1`
"""
function BesselA(
  m::Integer,
  x::Vector{<:Number},
  nu::Number
)
  if real(nu) <= -1
    throw(DomainError(nu, "The real part of `nu` is smaller than -1."))
  end
  p = length(x)
  b = nu + (p+1) / 2
  return hypergeomPQ(m, Int64[], [b], -x, 2) / mvgamma(b, p)
end

using ..HypergeomPQ
using ..Mvgamma

end # end module Bessel

using .HypergeomPQ
using .Mvgamma

end

ext: jl


text:
module Polynomial

import Base: print, println,  *, /, -, +, ==, divrem, div, rem, copy, one, zero, gcd, gcdx, lcm, call

using ResidueClass

export Poly, deg, ZeroPoly, Monom, lc, shift, Ring, EEA

## Comments ##
#
# No implicit type promotion to mymic closedness of mathematical operations in rings and fields
#
# Algorithms mostly based on Modern Computer Algebra (MCA) from Joachim von zu Gathen Jürgen Gerhard (Edition 2013)
#

immutable Poly{T}

    a::Vector{T}
    var::Symbol

    function Poly(a::Vector{T}, var::Symbol)
    	# delete leading zeros
    	i = length(a)
    	while a[i] == 0 && i > 1
    		i -= 1
    	end
    	i < length(a) ? new(a[1:i], var) : new(a,var)
    end
end

## Factory ##

Poly{T}(a::Vector{T}, var) = Poly{T}(a, symbol(var))
Poly{T}(a::Vector{T}) = Poly{T}(a, :X)

# returns poly X^d
function Monom{T}(sample::T, d::Int, var)
	a = Array{T}(d+1)
	for k in 1:d
		a[k] = zero(sample)
	end
	a[d+1] = one(sample)
	return Poly(a, symbol(var))
end

copy(p::Poly) = Poly(p.a, p.var)
one{T}(p::Poly{T}) = Poly([one(p[0])], p.var)
zero{T}(p::Poly{T}) = Poly([zero(p[0])], p.var)

one{T}(t::Type{Poly{T}}) = Poly([one(T)], :x)
# zero{T}(t::Type{Poly{T}}) = ZeroPoly(T, :x)

## Evaluation ##

# with this can poly p be evaluated with p(arg)
function call{T}(p::Poly{T},arg::T)
	# implementation of Horners Rule
	n = length(p.a)
	val = p.a[n]
	for i in n-1:-1:1
    	val *= arg
    	val += p.a[i]
	end
	return val
end

## Wieso klappt das nicht?!?

# with this can poly p be evaluated with p(arg)
# function call{T}(p::Poly{ResidueClass.Residue{T}},arg::T)
# 	# implementation of Horners Rule
# 	n = length(p.a)
# 	val = p.a[n]
# 	for i in n-1:-1:1
#     	val *= arg
#     	val += p.a[i]
# 	end
# 	return val
# end

## Printing ##

function printExponent(io::IO, var, i::Number)

	if i==0
	elseif i == 1
		print(io, var)
	else
		print(io, var,"^",i)
	end
end

function print{T<:Number}(io::IO, p::Poly{T})

	n = length(p.a)
	printed = false
	_one = one(p[0])
	_zero = zero(p[0])
	for i in n:-1:1
		c = p.a[i]
		if i < n && c > _zero && printed
			print(io, "+")
		end
		if (c != _zero) && (c!= _one || i == 1)
			print(io,c)
			printed = true
		end
		if (c != _zero)
			printExponent(io, p.var, i-1)
			printed = true
		end
	end
	if !printed
		return  print(io, "0")
	end
end

function print{T}(io::IO, p::Poly{T})

	n = length(p.a)
	printed = false
	_one = one(p[0])
	_zero = zero(p[0])
	for i in n:-1:1
		c = p.a[i]
		if i < n && c != _zero && printed
			print(io, "+")
		end
		if (c != _zero) && (c != _one || i == 1)
			print(io,"(",c,")")
			printed = true
		end
		if (c != _zero )
			printExponent(io, p.var, i-1)
			printed = true
		end
	end
	if !printed
		return  print(io, "0")
	end
end

## Indexing ##

getindex{T}(p::Poly{T}, i) = (i+1 > length(p.a) ? zero(p[0]) : p.a[i+1])

function setindex!{T}(p::Poly{T}, x, i0::Real)

	delta = i0+1 - length(p.a)
	if delta > 0
		_zero = zero(p[0])
		extension = Array{T}(delta)
		for k in 1:delta
			extension[k] = _zero
		end
		extension[delta] = convert(T,x)
		p.a = [p.a; extension]
	else
		p.a[i0+1] = convert(T,x)
	end
	return p
end

deg(p::Poly) = length(p.a) - 1

lc(p::Poly) = last(p.a)
sign(p::Poly) = sign(lc)

## Equates ##
function ==(p::Poly, q::Poly)
	if deg(p) != deg(q)
		return false
	elseif deg(p) > 0 && (p.var != q.var)
		return false
	else
		for i in 0:deg(p)
			if p[i] != q[i]
				return false
			end
		end
	end
	return true
end

==(a::WeakRef,p::Polynomial.Poly{WeakRef}) = deg(p) > 0 ? false : a == p[0]
==(p::Polynomial.Poly{WeakRef},a::WeakRef) = ==(a,p)

## Arithmetics ##

# plus #

function +{T,S}(p::Poly{T}, q::Poly{S})

	if p.var != q.var
		error("Can't add polynomials with var $(p.var) and $(q.var)")
	end

	return Poly([p[i] + q[i] for i=0:max(deg(p), deg(q))], p.var)
end

function +{T}(p::Poly{T}, s::T)

	if deg(p) == 0
		return p[0] + s
	end

	q = copy(p)
	q[0] += s
	return q
end

+{T}(s::T, p::Poly{T}) = +(p,s)

# minus #

function -{T}(p::Poly{T}, q::Poly{T})

	if p.var != q.var
		error("Can't add polynomials with var $(p.var) and $(q.var)")
	end

	return Poly([p[i] - q[i] for i=0:max(deg(p), deg(q))], p.var)
end

function -(p::Poly, s::Number)

	if deg(p) == 0
		return p[0] - s
	end

	q = copy(p)
	q[0] -= s
	return q
end

-(s::Number, p::Poly) = -(p,s)

# unary minus #
-(p::Poly) = Poly(-p.a, p.var)

# multiply #

*{T}(a::T, p::Poly{T}) = Poly(a * p.a, p.var)
*{T}(p::Poly{T}, a::T) = *(a,p)

function *{T}(p::Poly{T}, q::Poly{T})

	if p.var != q.var
		error("Can't multiply polynomials with var $(p.var) and $(q.var)")
	end
	n = deg(p)
	m = deg(q)

	c = Array{T}(n+m+1)
	# because of residue class -> m has to be known
	for k=1:n+m+1
		c[k]=zero(p[0])
	end
	for i in 1:n+1
		if p[i-1] != 0
    		c[i:m+i] += p[i-1] * q.a
    	end
	end

	return Poly(c, p.var)

end

function ^{T}(p::Poly{T}, n::Int)

    if n < 0
    	error("Cannot raise an polynomial p to a negative power -n")
    elseif n == 0
    	return one(p)
    end
    b = copy(p)
    for k in digits(n,2)[2:end]
        b = (k == 1) ? (b * b) * p : b * b
    end
    return b
end

# division #

function divrem{T}(p::Poly{T}, q::Poly{T})

	n = deg(p)
	m = deg(q)

	m > n && return (ZeroPoly(T, q.var), q)

	q_coeff = Array{T}(n-m+1)
	for k in 1:n-m+1
		q_coeff[k] = zero(p[0])
	end
	r = copy(p)
	u = one(p[0]) / lc(q)
	try u = convert(T, u)
	catch
		error("For divrem(p,q) is lc(q)=$(lc(q)) not a unit in $(T)")
	end

	for i in n-m:-1:0
		if deg(r) == m+i
			q_coeff[i+1] = lc(r) * u
			r = r - shift(q_coeff[i+1] * q, i)
		else
            q_coeff[i+1] = zero(p[0])
		end
	end

	return (Poly(q_coeff, q.var), r)
end

div{T}(p::Poly{T}, q::Poly{T}) = divrem(p,q)[1]
rem{T}(p::Poly{T}, q::Poly{T}) = divrem(p,q)[2]

function /{T}(p::Poly{T}, q::Poly{T})
	(div, rem) = divrem(p,q)
	if deg(rem) > 0 || rem[0] == zero(p[0])
		error("$(q) is no divisor of $(p)!")
	else
		return div
	end
end

function /{T}(p::Poly{T}, a::T)
    new_p = copy(p)
    for i in 1:length(p.a)
        p.a[i] /= a
    end
    return p
end

# For positive Integers multiplication with X^i
# For negative Integers a multiplication with X^-i, whereby we drop coefficients with j-i < 0,
# e.g. i = -2: X^5 + X -> X^3
function shift{T}(p::Poly{T}, i::Int)

	if i > 0
		zeros = Array{T}(i)
		for k in 1:i
			zeros[k] = zero(p[0])
		end
		return Poly([zeros; p.a], p.var)
	else
		if abs(i) > deg(p)
			return ZeroPoly(p[0], p.var)
		end
		return Poly(p.a[i+1:length(p.a)], p.var)
	end
end

## help functions ##

function isMonom{T}(p::Poly{T})

    lc(p) != one(p[0]) && return false
	_zero = zero(p[0])
 	deg(p) == _zero && return true

	for i in deg(p)-1:-1:0
		p[i] != _zero && return false
	end

	return true
end

#######################################################################################################################
#	A L G O R I T H M S
#######################################################################################################################

#######################################################################################################################
# E X T E N D E D  E U C L I D E A N  A L G O R I T H M (EEA)  F O R  Z  A N D  k[X]
#
# Implementation after MCA Algoriithm 3.14
#
# Input: f,g in R = Z or k[X]
# Output: normalized gcd(f,g) and s,t in R, such that s * f + t * g = gcd(f,g) (Bezout - Coefficients)
#

lu{T}(a::Poly{T}) = (a == zero(a[0])) ? one(a[0]) : lc(a)
lu{T<:Int}(a::T) =  (a == 1) ? one(T) : sign(a)

normal(a::Int) = abs(a)
normal{T}(a::Poly{T}) = a / lu(a)

## EEA / gcdx in Julia
function gcdx{T}(f::Poly{T}, g::Poly{T})

	 _one = one(g[0])
    _zero = zero(f[0])

    rho0 = lu(f)
    rho1 = lu(g)
    r0 = normal(f)
    r1 = normal(g)
    s0 = (rho0 == _one || rho0 == -_one) ? rho0 : rho0^-1
    s1 = _zero
    t0 = _zero
    t1 = (rho1 == _one || rho1 == -_one) ? rho1 : rho1^-1

    while r1 != _zero
        q = div(r0, r1)
        println("q:",q)
        rho2 = lu(r0 - q * r1)
        if rho2 == _zero
        	break
        end
        # eventually gcd in Z[x] possible, if-else necessary since 2 / 1 = 2.0 -> InexactError
        if (rho0 == _one || rho0 == -_one)
        	r2 = (r0 - q*r1) * rho2
        	s2 = (s0 - q*s1) * rho2
        	t2 = (t0 - q*t1) * rho2
        	r2 = (r0 - q*r1) * rho2
        else
        	s2 = (s0 - q*s1) / rho2
        	t2 = (t0 - q*t1) / rho2
        	r0 = r1
        	r1 = r2
        end
        rho0 = rho1
        rho1 = rho2
        t0 = t1
        t1 = t2
        s0 = s1
        s1 = s2
    end

    return (r0, s0, t0)
end

# efficient implementation if gcdx final
gcd{T}(f::Poly{T}, g::Poly{T}) = gcdx(f,g)[1]

lcm{T}(f::Poly{T}, g::Poly{T}) = (f * g) / gcd(f,g)

###########
##
## Interpolation

## Newton Interpolation
#  Implementation after Zippel, Effective Polynomial Computing (1993), p.217
function NewtonInterpol{T}(p::Poly{T}, values::Vector{T})

	(deg(p) != length(values)-1) &&	error("Not enough points to interpolate")

	f = zero(p)
	q = one(p)
	for a in values
		q_a = q(a)
		f = q_a != 0 ? f + (q_a^-1) * q * (p(a)-f(a)) : f
		q = Poly([-a, one(a)], p.var) * q
	end
	return f
end

end

max_stars_repo_name: saschatimme/Polynomial

max_issues_repo_name: saschatimme/Polynomial

max_forks_repo_name: saschatimme/Polynomial


text:
#
# qubit
#

function ITensors.space(::SiteType"qubit";
                        conserve_qns::Bool = false)
  if conserve_qns
    return [QN() => 2]
  end
  return 2
end

ITensors.state(::SiteType"qubit",
               ::ITensors.StateName"0") = 1

ITensors.state(::SiteType"qubit",
               ::ITensors.StateName"1") = 2

ITensors.op(::OpName"Id",
            ::SiteType"qubit",
            s::Index) =
  itensor([1 0; 0 1], s', dag(s))

ITensors.op(::OpName"X",
            ::SiteType"qubit",
            s::Index) =
  itensor([0.0 1.0; 1.0 0.0], s', dag(s))

ITensors.op(::OpName"iY",
            ::SiteType"qubit",
            s::Index) =
  itensor([0 1; -1 0], s', dag(s))

ITensors.op(::OpName"Y",
            ::SiteType"qubit",
            s::Index) =
  itensor([0 -im; im 0], s', dag(s))

ITensors.op(::OpName"Z",
            ::SiteType"qubit",
            s::Index) =
  itensor([1 0; 0 -1], s', dag(s))

ITensors.op(::OpName"H",
            ::SiteType"qubit",
            s::Index) =
  itensor([1/sqrt(2) 1/sqrt(2);
           1/sqrt(2) -1/sqrt(2)], s', dag(s))

ITensors.op(::OpName"Rx",
            ::SiteType"qubit",
            s::Index; θ::Number) =
  itensor([cos(θ/2)      -im*sin(θ/2.);
           -im*sin(θ/2.) cos(θ/2.)], s', dag(s))

ITensors.op(::OpName"Sw",
            ::SiteType"qubit",
            s1::Index,
            s2::Index) =
  itensor([1 0 0 0;
           0 0 1 0;
           0 1 0 0;
           0 0 0 1],
          s1', s2',
          dag(s1), dag(s2))

ITensors.op(::OpName"rand",
            ::SiteType"qubit",
            s::Index...) =
  randomITensor(prime.(s)...,
                dag.(s)...)

ITensors.op(::OpName"Cx",
            ::SiteType"qubit",
            s1::Index,
            s2::Index) =
  itensor([1 0 0 0;
           0 1 0 0;
           0 0 0 1;
           0 0 1 0],
          s1', s2',
          dag(s1), dag(s2))

ITensors.op(::OpName"T",
            ::SiteType"qubit",
            s1::Index,
            s2::Index,
            s3::Index) =
  itensor([1 0 0 0 0 0 0 0;
           0 1 0 0 0 0 0 0;
           0 0 1 0 0 0 0 0;
           0 0 0 1 0 0 0 0;
           0 0 0 0 1 0 0 0;
           0 0 0 0 0 1 0 0;
           0 0 0 0 0 0 0 1;
           0 0 0 0 0 0 1 0],
          s1', s2', s3',
          dag(s1), dag(s2), dag(s3))

ITensors.op(::OpName"noise",
            ::SiteType"qubit",
            s::Index...;
            krausind = hasqns(s[1]) ? Index([QN() => 2], "kraus") : Index(2, "kraus")) =
  randomITensor(prime.(s)..., dag.(s)..., krausind)


text:
using BinaryBuilder

name = "lapackpp"
version = v"2021.04.00"

# Collection of sources required to build PETSc. Avoid using the git repository, it will
# require building SOWING which fails in all non-linux platforms.
sources = [
    GitSource("https://bitbucket.org/icl/lapackpp.git", "31d969200a9f65390f56ac2ea48888bd10a13397")
]

# Bash recipe for building across all platforms

# Needs to add -Dcapi eventually once it's added to the cmake build system. Note yet available under CMAKAE toolchain.
script = raw"""
cd lapackpp
mkdir build && cd build
cmake \
  -DCMAKE_INSTALL_PREFIX=${prefix} \
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
  -DCMAKE_BUILD_TYPE="Release" \
  -Drun_result="0" \
  -Drun_result__TRYRUN_OUTPUT="ok" \
  -Dbuild_tests=no \
  -Duse_cmake_find_lapack=yes \
  ..
make -j${nproc}
make install
"""

# We attempt to build for all defined platforms
platforms = expand_cxxstring_abis(supported_platforms())
products = [
    LibraryProduct("liblapackpp", :liblapackpp)
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS32_jll"),
    Dependency("blaspp_jll"; compat="2021.04.02")
]
# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"7")

max_stars_repo_name: jakebolewski/Yggdrasil

max_issues_repo_name: jakebolewski/Yggdrasil

max_forks_repo_name: jakebolewski/Yggdrasil


text:
"""
    loglik_obs!(d, y, μ, wt, ϕ)
Get the loglikelihood from the GLM.jl package for each observation
"""
function loglik_obs end

loglik_obs(::Bernoulli, y, μ, wt, ϕ) = wt*GLM.logpdf(Bernoulli(μ), y)
loglik_obs(::Binomial, y, μ, wt, ϕ) = GLM.logpdf(Binomial(Int(wt), μ), Int(y*wt))
loglik_obs(::Gamma, y, μ, wt, ϕ) = wt*GLM.logpdf(Gamma(inv(ϕ), μ*ϕ), y)
loglik_obs(::InverseGaussian, y, μ, wt, ϕ) = wt*GLM.logpdf(InverseGaussian(μ, inv(ϕ)), y)
loglik_obs(::Normal, y, μ, wt, ϕ) = wt*GLM.logpdf(Normal(μ, sqrt(abs(ϕ))), y)
loglik_obs(::Poisson, y, μ, wt, ϕ) = logpdf(Poisson(μ), y)

# this gets the loglikelihood from the glm.jl package for the component density
"""
    component_loglikelihood!(gc::GLMCopulaVCObs{T, D, Link})
Calculates the loglikelihood of observing `y` given mean `μ`, a distribution
`d` using the GLM.jl package.
"""
function component_loglikelihood(gc::Union{GLMCopulaVCObs{T, D, Link}, GLMCopulaARObs{T, D, Link}, GLMCopulaCSObs{T, D, Link}}) where {T <: BlasReal, D, Link}
  logl = zero(T)
    @inbounds for j in 1:gc.n
      logl += QuasiCopula.loglik_obs(gc.d, gc.y[j], gc.μ[j], gc.wt[j], 1.0)
  end
  logl
end

### loglikelihood functions
"""
    component_loglikelihood!(gc::Union{NBCopulaVCObs{T, D, Link}, NBCopulaARObs{T, D, Link}, NBCopulaCSObs{T, D, Link}}, r::T)
Calculates the loglikelihood of observing `y` given parameters for `μ` and `r` for Negative Binomial distribution using the GLM.jl package.
"""
function component_loglikelihood(gc::Union{NBCopulaVCObs{T, D, Link}, NBCopulaARObs{T, D, Link}, NBCopulaCSObs{T, D, Link}}, r::T) where {T <: BlasReal, D<:NegativeBinomial{T}, Link}
    logl = zero(T)
    @inbounds for j in 1:gc.n
        logl += logpdf(D(r, r / (gc.μ[j] + r)), gc.y[j])
    end
    logl
end

"""
    component_loglikelihood!(gc::Poisson_Bernoulli_VCObs{T, VD, VL})
Calculates the loglikelihood of observing `y` given mean `μ`, a distribution
`d` with mixed types of poisson and bernoulli distributions using the GLM.jl package.
"""
function component_loglikelihood(gc::Poisson_Bernoulli_VCObs{T, VD, VL}) where {T <: BlasReal, VD, VL}
  logl = zero(T)
    @inbounds for j in 1:gc.n
      logl += QuasiCopula.loglik_obs(gc.vecd[j], gc.y[j], gc.μ[j], gc.wt[j], 1.0)
  end
  logl
end


text:
# Interface -------

function forks(owner, repo; auth = AnonymousAuth(), options...)
    forks(auth, owner, repo; options...)
end

function forks(auth, owner, repo; headers = Dict(), options...)
    authenticate_headers(headers, auth)
    r = get(URI(API_ENDPOINT; path = "/repos/$owner/$repo/forks");
            headers = headers,
            options...)

    handle_error(r)

    data = JSON.parse(r.data)
    map!( r -> Repo(r), data)
end

function fork(owner, repo, organization = ""; auth = AnonymousAuth(), options...)
    fork(auth, owner, repo, organization; options...)
end

function fork(auth, owner, repo, organization = ""; headers = Dict(),
                                                    json = Dict(),
                                                    options...)
    authenticate_headers(headers, auth)

    if organization != ""
        json["organization"] = organization
    end

    r = post(URI(API_ENDPOINT; path = "/repos/$owner/$repo/forks");
            headers = headers,
            json = json,
            options...)

    handle_error(r)

    data = JSON.parse(r.data)
end


text: println("hello, " * "world")


text:
function init_pseudo_observable_mappings!(m::Model1002)

    pseudo_names = [:y_t, :y_f_t, :NaturalRate, :π_t, :OutputGap, :ExAnteRealRate, :LongRunInflation,
                    :MarginalCost, :Wages, :FlexibleWages, :Hours, :FlexibleHours, :z_t,
                    :Expected10YearRateGap, :NominalFFR, :Expected10YearRate,
                    :Expected10YearNaturalRate,
                    :ExpectedNominalNaturalRate, :NominalRateGap, :LaborProductivityGrowth,
                    :u_t, :NominalWageGrowth]

    if subspec(m) == "ss12"
        to_add = [:g_t, :b_t, :μ_t, :λ_f_t, :λ_w_t, :rm_t, :σ_ω_t, :μ_e_t,
                  :γ_t, :π_star_t, :e_lr_t, :e_tfp_t, :e_gdpdef_t, :e_corepce_t, :e_gdp_t, :e_gdi_t]
        pseudo_names = vcat(pseudo_names, to_add)
    end

    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        push!(pseudo_names, :Sinf_t, :Sinf_w_coef_t, :ι_p, :πtil_t, :πtil_t1, :e_tfp_t)
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            push!(pseudo_names, :e_tfp_t1)
        end
    end

    if haskey(m.settings, :add_laborshare_measurement)
        if get_setting(m, :add_laborshare_measurement)
            push!(pseudo_names, :laborshare_t)
        end
    end

    # Create PseudoObservable objects
    pseudo = OrderedDict{Symbol,PseudoObservable}()
    for k in pseudo_names
        pseudo[k] = PseudoObservable(k)
    end

    # Fill in names and reverse transforms
    pseudo[:y_t].name = "Output Growth"
    pseudo[:y_t].longname = "Output Growth Per Capita"

    pseudo[:y_f_t].name = "Flexible Output Growth"
    pseudo[:y_f_t].longname = "Output that would prevail in a flexible-price economy."

    pseudo[:NaturalRate].name = "Real Natural Rate"
    pseudo[:NaturalRate].longname = "The real interest rate that would prevail in a flexible-price economy."
    pseudo[:NaturalRate].rev_transform = quartertoannual

    pseudo[:π_t].name = "Inflation"
    pseudo[:π_t].longname = "Inflation"
    pseudo[:π_t].rev_transform = quartertoannual

    pseudo[:OutputGap].name = "Output Gap"
    pseudo[:OutputGap].longname = "Output Gap"

    pseudo[:ExAnteRealRate].name = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].longname = "Ex Ante Real Rate"
    pseudo[:ExAnteRealRate].rev_transform = quartertoannual

    pseudo[:LongRunInflation].name = "Long Run Inflation"
    pseudo[:LongRunInflation].longname = "Long Run Inflation"
    pseudo[:LongRunInflation].rev_transform = quartertoannual

    pseudo[:MarginalCost].name = "Marginal Cost"
    pseudo[:MarginalCost].longname = "Marginal Cost"

    pseudo[:Wages].name = "Wages"
    pseudo[:Wages].longname = "Wages"

    pseudo[:FlexibleWages].name = "Flexible Wages"
    pseudo[:FlexibleWages].longname = "Wages that would prevail in a flexible-wage economy"

    pseudo[:Hours].name = "Hours"
    pseudo[:Hours].longname = "Hours"

    pseudo[:FlexibleHours].name     = "Flexible Hours"
    pseudo[:FlexibleHours].longname = "Flexible Hours"

    pseudo[:z_t].name     = "z_t (Technology Growth minus Steady State Growth)"
    pseudo[:z_t].longname = "z_t (Technology Growth minus Steady State Growth)"

    pseudo[:Expected10YearRateGap].name     = "Expected 10-Year Rate Gap"
    pseudo[:Expected10YearRateGap].longname = "Expected 10-Year Rate Gap"
    pseudo[:Expected10YearRateGap].rev_transform = quartertoannual

    pseudo[:NominalFFR].name     = "Nominal FFR"
    pseudo[:NominalFFR].longname = "Nominal FFR at an annual rate"
    pseudo[:NominalFFR].rev_transform = quartertoannual

    pseudo[:Expected10YearRate].name     = "Expected 10-Year Rate"
    pseudo[:Expected10YearRate].longname = "Expected 10-Year Interest Rate"
    pseudo[:Expected10YearRate].rev_transform = quartertoannual

    pseudo[:Expected10YearNaturalRate].name     = "Expected 10-Year Natural Rate"
    pseudo[:Expected10YearNaturalRate].longname = "Expected 10-Year Natural Rate of Interest"
    pseudo[:Expected10YearNaturalRate].rev_transform = quartertoannual

    pseudo[:ExpectedNominalNaturalRate].name     = "Expected Nominal Natural Rate"
    pseudo[:ExpectedNominalNaturalRate].longname = "Natural Rate + Expected Inflation"
    pseudo[:ExpectedNominalNaturalRate].rev_transform = quartertoannual

    pseudo[:NominalRateGap].name     = "Nominal Rate Gap"
    pseudo[:NominalRateGap].longname = "Nominal FFR - Nominal Natural Rate"
    pseudo[:NominalRateGap].rev_transform = quartertoannual

    pseudo[:LaborProductivityGrowth].name     = "Labor Productivity Growth"
    pseudo[:LaborProductivityGrowth].longname = "Labor Productivity Growth Rate"
    pseudo[:LaborProductivityGrowth].rev_transform = quartertoannual

    pseudo[:u_t].name     = "u_t"
    pseudo[:u_t].longname = "u_t"

    pseudo[:NominalWageGrowth].name = "Nominal Wage Growth"
    pseudo[:NominalWageGrowth].longname = "Nominal Wage Growth"

    if subspec(m) in ["ss13", "ss14", "ss15", "ss16", "ss17", "ss18", "ss19"]
        pseudo[:Sinf_t].name     = "Sinf_t"
        pseudo[:Sinf_t].longname = "Sinf_t, PDV of Emc_t"
        pseudo[:Sinf_w_coef_t].name     = "Sinf_w_coef_t"
        pseudo[:Sinf_w_coef_t].longname = "Sinf_w_coef_t, PDV of Emc_t multiplied by coefficient"
        pseudo[:ι_p].name = "iota_p"
        pseudo[:ι_p].longname = "iota_p"
        pseudo[:πtil_t].name     = "pitil_t"
        pseudo[:πtil_t].longname = "Fundamental Inflation"
        pseudo[:πtil_t1].name     = "pitil_t1"
        pseudo[:πtil_t1].longname = "Fundamental Inflation Lag 1"
        pseudo[:e_tfp_t].name = "e_tfp_t"
        pseudo[:e_tfp_t].longname = "e_tfp_t"
        if subspec(m) in ["ss14", "ss15", "ss16", "ss18", "ss19"]
            pseudo[:e_tfp_t1].name = "e_tfp_t1"
            pseudo[:e_tfp_t1].longname = "e_tfp_t1"
        end
    end

    # Other exogenous processes
    if subspec(m) == "ss12"
        for i in to_add
            pseudo[i].name = string(DSGE.detexify(i))
            pseudo[i].longname = string(DSGE.detexify(i))
        end
    end

    if haskey(m.settings, :add_laborshare_measurement)
        if get_setting(m, :add_laborshare_measurement)
            pseudo[:laborshare_t].name     = "Log Labor Share"
            pseudo[:laborshare_t].longname = "Log Labor Share"
        end
    end

    # Add to model object
    m.pseudo_observable_mappings = pseudo
end


text:
using Statistics
using StatsBase
using SpecialFunctions
using Roots

include("load_data.jl")

###############################################################################
###############################################################################
### functions that work on interaction activity

"""
    coefficient_variation(cas::contact_activities{T}; jackknife_error=false) where T

calculates the coefficient of variation as the ratio of std deviation to the
mean (sigma/mu) accross all individuals from the contact_activities object.
This means that the mean and the variance are really calculates over all
individuals.

With the option `jackknife_error=true` a jackknife error is computed and
returned (which is factor 100+ faster than our cas convenience solution because
we can eaily precompute the partial contributions for the jackknife procedure)
"""
function coefficient_variation(cas::contact_activities{T}; jackknife_error=false) where T
    # data = vector of tuples (mean(a), mean(a^2)) for local activities a per id
    # avg = 1/N sum(mean(a))
    # var = 1/N (sum(mean(a^2)) - avg^2
    function estimate_cv(data)
        N = length(data)
        avg = mean(getindex.(data,1))
        std = sqrt((sum(getindex.(data,2)) - N*avg^2)/(N-1))
        return std/avg
    end
    #print("... coefficient of variation:\n") -> move to analysis all functions
    data = Vector{Tuple{Float64, Float64}}(undef, length(activities(cas)))
    @showprogress 1 for (i,a) in enumerate(activities(cas))
        data[i] = (mean(a), mean(x->x^2, a))
    end
    cv_est = estimate_cv(data)
    if jackknife_error
        return jackknife(estimate_cv, data, naive=cv_est)
    else
        return cv_est
    end
end

"""
    autocorrelation_function(cas::contact_activities, lags;)

calculate the autocorrelation function of the interaction activities accross
participants for the provided lags.
"""
function autocorrelation_function(cas::contact_activities{I,T}, lags::Union{AbstractRange{Int64}, AbstractVector{Int64}}) where {I,T}
    @assert length(lags) > 0
    avg = 0.0
    for a in activities(cas)
        avg += sum(a)/length(a)
    end
    avg /= length(activities(cas))
    #avg = sum(sum.(activities(cas)))/length(activities(cas))/length(times(cas))
    cov = zeros(Float64, length(lags))
    p = Progress(length(cas), 1, "Autocorrelation: ", offset=1)
    for a in activities(cas)
        for (index, k) in enumerate(lags)
            for i in 1:length(a)-k
                cov[index] += a[i]*a[i+k]
            end
        end
        next!(p)
    end
    cov ./= length(activities(cas)).*(length(times(cas)).-lags)

    C  = cov .- avg^2
    C ./= C[1]
    return C
end

"""
    rate(cas::contact_activities, time_range)

calculate the rate of interactions for the times in `time_range`
(start:step:end) averaged over individuals and repetitions of `time_range` in
the data (e.g. if time_range is 0:timestep:1week) and the data has 4
wees, then the average ist also across weeks.
"""
function rate(cas::contact_activities{I,T}, time_range::AbstractRange) where {I,T}
    @assert step(time_range) >= step(times(cas))

    rate = Histogram(time_range)
    norm = Histogram(time_range)
    time_max = last(time_range)
    for a in activities(cas)
        for (i, t) in enumerate(times(cas))
            # modulo for periodic boundaries
            rate[ t%time_max ] += a[i]
            norm[ t%time_max ] += 1
        end
    end

    # rate is avg. number of interactions per timestep
    rate = float(rate)
    rate.weights ./= norm.weights

    return rate
end

###############################################################################
###############################################################################
### functions that operate on list of durations
"""
returns the distribution of contact durations from `list_durations`, which is as usual list of lists of durations

# Example for data preparation:
```
    list_durations = [getindex.(encounter, 4) for encounter in list_encounter];
```

# Error
If interested in the statistics of contacts (most natural) one has to estimate the error of the distribution weights.
Important is that the jackknife routine does not comply with leaving out m=0 elements
```
    full_dist = distribution_contact_durations(list_durations, timestep=timestep(experiment));
    P, Pj, Perr = jackknife_mj(x->distribution_contact_durations(x, edges=full_dist.edges[1]).weights, list_durations[length.(list_durations) .> 0], naive=full_dist.weights);
```

In case one is interested in logartihmic binning, this should be the way to achieve it
```
    xbin, Pbin = logbin(full_dist)
    xP, xP_J, xP_err = jackknife_mj(x->hcat(logbin(distribution_contact_durations(x, edges=full_dist.edges[1]))...)[:,1:2], list_durations[length.(list_durations) .> 0], naive=hcat(xbin,Pbin));
```
which equals
```
    xbin_n, xbin_J, xbin_err = jackknife_mj(x->logbin(distribution_contact_durations(x, edges=full_dist.edges[1]))[1], list_durations[length.(list_durations) .> 0], naive=xbin);
    Pbin_n, Pbin_J, Pbin_err = jackknife_mj(x->logbin(distribution_contact_durations(x, edges=full_dist.edges[1]))[2], list_durations[length.(list_durations) .> 0], naive=Pbin);
```

"""
function distribution_durations(list_x::Vector{T}; edges=missing, timestep=missing) where T
    # flatten the list of lists (not very expansive as it turns out)
    flat_list = vcat(list_x...)
    # construct edges of distribution if not provided
    if ismissing(edges)
        if ismissing(timestep)
            throw(error("either distribution `edges` or `timestep` of experiment have to be specified via keyword arguments"))
        else
            edges = 0:timestep:maximum(flat_list)+timestep
        end
    end
    # sort elements of list into Histogram
    dist = fit(Histogram{Float64}, flat_list, edges)
    normalize!(dist)
    return dist
end

###############################################################################
###############################################################################
### functions that operate on encounter trains (or derivatives such as list_dts)
"""
    rate(ets::encounter_trains, time_range)

calculate the rate of contacts for the times in `time_range`
(start:step:end) averaged over individuals and repetitions of `time_range` in
the data (e.g. if time_range is 0:timestep:1week) and the data has 4
weeks, then the average ist also across weeks.

# Error estimation
should be based on ordinary jackknife method, because the quantity of rate has
same number of entries for each train
 ```
    r,rj,err = jackknife(x->rate(x, 0:300:seconds_from_days(7)), ets)
```

"""
function rate(
        ets::encounter_trains{I,T},
        time_range::AbstractRange;
        # optional arguments
    ) where {I,T}
    @assert step(time_range) >= timestep(ets)
    # for now because I am unsure how do deal with periodic ranges that do not
    # start at 0 -> easier to align data acordingly
    @assert first(time_range) == 0

    num_contacts = Histogram(time_range)
    time_max = last(time_range)
    for train in trains(ets)
        for t in train
            # modulo for periodic boundaries
            num_contacts[ t%time_max ] += 1
        end
    end

    # normalization considers how good time_range fits in duration(ets)
    norm = Histogram(time_range)
    time_max = last(time_range)
    norm.weights .+= floor(Int, duration(ets) / time_max) * length(trains(ets))
    time_left = duration(ets) % time_max
    if time_left > 0
        norm.weights[1:StatsBase.binindex(norm, time_left)] .+= length(trains(ets))
    end

    # rate is avg. number of contacts per second (hence need to devide
    # additionally by time step)
    rate = float(num_contacts)
    rate.weights ./= norm.weights
    rate.weights ./= step(time_range)

    return rate
end

"""

returns the rate of encounters in `time_range` (with according timesteps, unit
of dataset) conditioned on having a contact at time 0.
"""
function conditional_encounter_rate(
        ets::encounter_trains{I,T},
        time_range::AbstractRange;
        # optional arguments
    ) where {I,T}
    num_contacts = Histogram(time_range)
    norm = Histogram(time_range)

    for train in trains(ets)
        #collect statistic over all contacts within the train
        for (i, contact_time) in enumerate(train)
            # add next contacts to histogram
            for j in i+1:length(train)
                time = @inbounds train[j]
                push!(num_contacts, time - contact_time)
                #if time - contact_time > last(num_contacts.edges[1])
                #    break
                #end
            end

            # add a 1 to all those bins that are covered by the dataset
            time_max = duration(ets) - contact_time
            #binindex returns max bin if time_max crosses boundaries
            index_max = StatsBase.binindex(norm, time_max)
            if index_max > length(norm.weights)
                @inbounds norm.weights .+= 1
            else
                for j in 1:index_max
                    @inbounds norm.weights[j] += 1
                end
            end
        end
    end

    rate = float(num_contacts)
    rate.weights ./= step(time_range)
    rate.weights ./= norm.weights

    return rate
end

function integrate(rate::AbstractHistogram, interval::Tuple)
    @assert(first(interval) <= last(interval))
    @assert(first(interval) >= first(rate.edges[1]))
    @assert(first(interval) < last(rate.edges[1]))

    first_idx = StatsBase.binindex(rate, first(interval))
    last_idx = StatsBase.binindex(rate, last(interval))

    integral = sum(rate.weights[first_idx:last_idx])

    return integral*step(rate.edges[1])
end

"""
    distribution_number_encoutner(ets, window, [pattern,])

calculates the distribution of the number of contacts within a certain time
window. With the keyword `pattern` one may specify a pattern of reoccuring
time-windows that should be distinguished. Using pattern one should be aware of
the reference time of the experiment (which can be to some extent be controlled
with the keyword offset). (Copenhagen starts on Sunday assumed at 0:00)

Examples:
calculate distribution of daily number of contacts irrespective of which day
```
    distribution_number_encoutner(ets, seconds_from_days(1))
``

calculate distribution of daily number of contacts for each day of the week separately
```
    distribution_number_encoutner(ets, seconds_from_days(1), pattern=1:7)
``

calculate distribution of daily number of contacts distinguishing weekdays and weekends
```
    distribution_number_encoutner(ets, seconds_from_days(1), pattern=[1,2,2,2,2,2,1])
``

"""
function distribution_number_encounter(
        ets::encounter_trains{I,T},
        window::Real,
        pattern::Union{AbstractRange, AbstractVector};
        ### optional
        offset::Real = 0,
        edges=missing,
    ) where {I,T}
    # find unique elements of pattern
    unique_labels = sort!(unique(pattern))

    # create Histogram to sort count contacts
    num_contacts = Histogram(0:window:duration(ets))

    # create dictionary to sort the number of contacts
    list_num = Dict(unique_labels .=> [Int64[] for i in 1:length(unique_labels)])

    # iterate over contact trains and sort into lists (each part of pattern is
    # `window` long)
    p = Progress(length(ets), 1, "Distribution: ", offset=1)
    for train in trains(ets)
        next!(p)
        num_contacts.weights .= 0
        for time in train
            push!(num_contacts, time-offset)
        end
        # add to list
        for (i, num) in enumerate(num_contacts.weights)
            push!(list_num[pattern[1+(i-1)%length(pattern)]], num)
        end
    end

    if ismissing(edges)
        list_dist = [fit(Histogram{Float64}, list_num[label], 0:1:maximum(list_num[label])+1) for label in unique_labels]
    else
        list_dist = [fit(Histogram{Float64}, list_num[label], edges) for label in unique_labels]
    end
    normalize!.(list_dist)

    return Dict(unique_labels .=> list_dist)
end
distribution_number_encounter(ets::encounter_trains{I,T}, window::Real; offset::Real = 0, edges=missing, return_type=Dict) where {I,T} = distribution_number_encounter(ets, window, [1,], offset=offset, edges=edges)

"""
    inter_encounter_intervals(ets)

calculate the list of inter-encoutner times for each id in `ets` (of type encoutner
trains) and return as vector of vector of inter-encoutner times
(Vector{Vector{T}}).
"""
function inter_encounter_intervals(ets::encounter_trains{I,T})::Vector{Vector{T}} where {I,T}
    list_dts = Vector{T}[]
    for train in trains(ets)
        dts = diff(train)
        push!(list_dts, dts)
    end
    return list_dts
end

function inter_encounter_intervals(trains::Vector{Vector{T}})::Vector{Vector{T}} where {I,T}
    list_dts = Vector{T}[]
    for train in trains
        dts = diff(train)
        push!(list_dts, dts)
    end
    return list_dts
end

"""
    autocorrelation_function(list_dts, lags)

calculate autocorrelation function of inter-contact intervals. Due to the
finite-sample bias, we precalculate the stationary mean interval for each lag
once for the x-data set and once for the y-data set, cf. [Spitzner et al., PLOS
One (2021)].

# Error estimation:
The statistics of interest here is the inter-contact intervals (dts), which
differe from list to list. Hence, jackknife_mj should be used:
```
    c, cj, cerr = jackknife_mj(x->autocorrelation_function(x, lags), list_dts)
```
"""
function autocorrelation_function(list_dts::Vector{Vector{T}}, lags::Union{AbstractRange{Int64}, AbstractVector{Int64}}) where T
    @assert length(lags) > 0

    # calculate global mean assuming stationary inter-contact intervals across
    # individuals
    sum_dts = zero(T)
    num_dts = zero(Int64)
    for dts in list_dts
        sum_dts += sum(dts)
        num_dts += length(dts)
    end
    mean_dt = float(sum_dts)/float(num_dts)
    #mean_dt = sum(sum.(list_dts))/sum(length.(list_dts))

    # calculate actual autocorrelation function E[ (X_1 - E[X_1])(X_2 - E[X_2]) ]
    C  = zeros(Float64, length(lags))
    N  = zeros(Int64, length(lags))
    p = Progress(length(list_dts), 1, "Autocorrelation: ", offset=1)
    for dts in list_dts
        for (index,k) in enumerate(lags)
            for i in 1:length(dts)-k
                term = (dts[i] - mean_dt) * (dts[i+k] - mean_dt)
                C[index] += term
            end
            # N is iteratively calculated because lengths(dts) is not the same
            # across elements
            N[index] += max(length(dts)-k, 0)
        end
        next!(p)
    end

    C  ./= N
    C  ./= C[1]
    return C
end
#wrapper functions
autocorrelation_function(dts::Vector{T}, lags) where T = autocorrelation_function([dts,], lags)
autocorrelation_function(ets::encounter_trains, lags) = autocorrelation_function(inter_contact_times(ets), lags)

function autocorrelation_function_hack(list_dts::Vector{Vector{T}}, lags) where T
    # precalculate mean dt for each lag across trains, distinguishing between
    # [1:end-k] and [1+k:end]
    mean_dt_x = zeros(length(lags))
    mean_dt_y = zeros(length(lags))
    N  = zeros(length(lags))
    for dts in list_dts
        for (j,k) in enumerate(lags)
            mean_dt_x[j] += sum(dts[1:end-k])
            mean_dt_y[j] += sum(dts[1+k:end])
            N[j] += max(length(dts)-k, 0)
        end
    end
    mean_dt_x ./= N
    mean_dt_y ./= N

    # calculate actual autocorrelation function E[ (X_1 - E[X_1])(X_2 - E[X_2]) ]
    C  = zeros(length(lags))
    C2 = zeros(length(lags))
    N  = zeros(length(lags))
    for dts in list_dts
        for (j,k) in enumerate(lags)
            for i in 1:length(dts)-k
                term = (dts[i] - mean_dt_x[j])*(dts[i+k] - mean_dt_y[j])
                C[j]  += term
            end
            N[j] += max(length(dts)-k, 0)
        end
    end
    C  ./= N
    Var = C[1]
    C  ./= Var

    return C
end

###############################################################################
###############################################################################
### cluster

"""
sample encounter cluster from the list of inter-encounter intervals
and a minimum separation that specifies that larger inter-encoutner intervals terminate the cluster

minimum_separation should be on the order of mean(duration)~8*300 for CNS
"""
function sample_cluster_sizes(list_dts::Vector{Vector{T}}, minimum_separation::Number) where T
    samples = [Int[] for i in 1:length(list_dts)]
    for (i,dts) in enumerate(list_dts)
        # start from cluster size 1 because we work on level of dts
        size = 1
        for dt in dts
            if dt < minimum_separation
                size += 1
            else
                push!(samples[i], size)
                size = 1
            end
        end
    end

    return samples
end

"""
evaluate distribution of contact cluster from contact activities as the sum of
successive simultaneious contacts.
"""
function sample_cluster_sizes(cas::contact_activities{I,T}) where {I,T}
    samples = [Int[] for i in 1:length(cas)]
    for (i,activity) in enumerate(activities(cas))
        size = 0
        for a in activity
            if size > 0 && a==0
                push!(samples[i], size)
                size = 0
            end
            size += a
        end
    end

    return samples
end

"""
    This could also be used for distribution_duration
"""
function distribution(samples::Vector{T}; edges=missing, resolution=1) where T
    # flatten the list of lists (not very expansive as it turns out)
    flat_samples = vcat(samples...)
    if ismissing(edges)
        dist = fit(Histogram{Float64}, flat_samples, 0:resolution:maximum(flat_samples)+resolution)
    else
        dist = fit(Histogram{Float64}, flat_samples, edges)
    end
    normalize!(dist)

    return dist
end


text:
module KroneckerTools

using LinearAlgebra
using LinearAlgebra.BLAS
using QuasiTriangular
import QuasiTriangular: A_mul_B!, At_mul_B!, A_mul_Bt!
import Base.convert

include("ExtendedMul.jl")
export a_mul_kron_b!, a_mul_b_kron_c!, kron_at_kron_b_mul_c!, a_mul_b_kron_c_d!, at_mul_b_kron_c!, a_mul_b_kron_ct!

"""
Content:
c = a*(b ⊗ b ⊗ ... ⊗ b)
a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, order::Int64)
a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::AbstractVector,work::AbstractVector)
d = a*b*(c ⊗ c ⊗ ... ⊗ c)
a_mul_b_kron_c!(d::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, order::Int64,work1::AbstractVector,work2::AbstractVector)
d = a^T*b*(c ⊗ c ⊗ ... ⊗ c)
at_mul_b_kron_c!(d::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
d = a*b*(c^T ⊗ c^T ⊗ ... ⊗ c^T)
a_mul_b_kron_ct!(d::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
e = a*b*(c ⊗ d ⊗ ... ⊗ d)
a_mul_b_kron_c_d!(e::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, d::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
d = (a^T ⊗ a^T ⊗ ... ⊗ a^T ⊗ b)*c
kron_at_kron_b_mul_c!(d::AbstractVector, offset_d::Int64, a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, offset_c::Int64, work1::AbstractVector, work2::AbstractVector, offset_w::Int64)
kron_at_kron_b_mul_c!(d::AbstractVector, a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, work1::AbstractVector, work2::AbstractVector)
kron_at_kron_b_mul_c!(a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, offset_c::Int64, work::AbstractVector)
kron_at_kron_b_mul_c!(a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, work::AbstractVector)
c = (a ⊗ a ⊗ ... ⊗ a)*b
kron_a_mul_b!(c::AbstractVector, a::AbstractMatrix, order::Int64, b::AbstractVector, q::Int64, work1::AbstractVector, work2::AbstractVector)
c = (I_p ⊗ a ⊗ I_q)*b
kron_mul_elem!(c::AbstractVector, offset_c::Int64, a::AbstractMatrix, b::AbstractVector, offset_b::Int64, p::Int64, q::Int64)
kron_mul_elem!(c::AbstractVector, a::AbstractMatrix, b::AbstractVector, p::Int64, q::Int64)
c = (I_p ⊗ a^T ⊗ I_q)*b
kron_mul_elem_t!(c::AbstractVector, offset_c::Int64, a::AbstractMatrix, b::AbstractVector, offset_b::Int64, p::Int64, q::Int64)
kron_mul_elem_t!(c::AbstractVector, a::AbstractMatrix, b::AbstractVector, p::Int64, q::Int64)
"""

"""
    kron_mul_elem!(p::Int64, q::Int64, m::Int64, a::AbstractMatrix, b::AbstractVector, c::AbstractVector)

Performs (I_p ⊗ a ⊗ I_q) b, where m,n = size(a). The result is stored in c.
"""
function kron_mul_elem!(c::AbstractVector, offset_c::Int64, a::AbstractMatrix, b::AbstractVector, offset_b::Int64, p::Int64, q::Int64)
    m, n = size(a)
    length(b) >= n*p*q || throw(DimensionMismatch("The dimension of vector b, $(length(b)) doesn't correspond to order, ($p, $q)  and the dimensions of the matrix, $(size(a))"))
    length(c) >= m*p*q || throw(DimensionMismatch("The dimension of the vector c, $(length(c)) doesn't correspond to order, ($p, $q)  and the dimensions of the matrix, $(size(a))"))

    begin
        if p == 1 && q == 1
            # a*b
            A_mul_B!(c, offset_c, a, b, offset_b, 1)
        elseif q == 1
            #  (I_p ⊗ a)*b = vec(a*[b_1 b_2 ... b_p])
            A_mul_B!(c, offset_c, a, b, offset_b, p)
        elseif p == 1
            # (a ⊗ I_q)*b = (b'*(a' ⊗ I_q))' = vec(reshape(b,q,m)*a')
            A_mul_Bt!(c, offset_c, b, offset_b, q, n, a)
        else
            # (I_p ⊗ a ⊗ I_q)*b = vec([(a ⊗ I_q)*b_1 (a ⊗ I_q)*b_2 ... (a ⊗ I_q)*b_p])
            mq = m*q
            nq = n*q
            for i=1:p
                A_mul_Bt!(c, offset_c, b, offset_b, q, n, a)
                offset_b += nq
                offset_c += mq
            end
        end
    end
end

function kron_mul_elem!(c::AbstractVector, a::AbstractMatrix, b::AbstractVector, p::Int64, q::Int64)
    kron_mul_elem!(c, 1, a, b, 1, p, q)
end

"""
    kron_mul_elem_t!(p::Int64, q::Int64, m::Int64, a::AbstractMatrix, b::AbstractVector, c::AbstractVector)

Performs (I_p ⊗ a' ⊗ I_q) b, where m,n = size(a). The result is stored in c.
"""
function kron_mul_elem_t!(c::AbstractVector, offset_c::Int64, a::AbstractMatrix, b::AbstractVector, offset_b::Int64, p::Int64, q::Int64)
    m, n = size(a)
    length(b) >= m*p*q || throw(DimensionMismatch("The dimension of vector b, $(length(b)) doesn't correspond to order, ($p, $q)  and the dimensions of the matrix, $(size(a))"))
    length(c) >= n*p*q || throw(DimensionMismatch("The dimension of the vector c, $(length(c)) doesn't correspond to order, ($p, $q)  and the dimensions of the matrix, $(size(a))"))

    begin
        if p == 1 && q == 1
            # a'*b
            At_mul_B!(c, offset_c, a, b, offset_b, 1)
        elseif q == 1
            #  (I_p ⊗ a')*b = vec(a'*[b_1 b_2 ... b_p])
            At_mul_B!(c, offset_c, a, b, offset_b, p)
        elseif p == 1
            # (a' ⊗ I_q)*b = (b'*(a ⊗ I_q))' = vec(reshape(b,q,m)*a)
            A_mul_B!(c, offset_c, b, offset_b, q, m, a)
        else
            # (I_p ⊗ a' ⊗ I_q)*b = vec([(a' ⊗ I_q)*b_1 (a' ⊗ I_q)*b_2 ... (a' ⊗ I_q)*b_p])
            mq = m*q
            nq = n*q
            for i=1:p
                A_mul_B!(c, offset_c, b, offset_b, q, m, a)
                offset_b += mq
                offset_c += nq
            end
        end
    end
end

function kron_mul_elem_t!(c::AbstractVector, a::AbstractMatrix, b::AbstractVector, p::Int64, q::Int64)
    kron_mul_elem_t!(c, 1, a, b, 1, p, q)
end

"""
    a_mul_kron_b!(c::AbstractVector, a::AbstractVecOrMat, b::AbstractMatrix, order::Int64)

Performs a*(b ⊗ b ⊗ ... ⊗ b). The solution is returned in matrix c. order indicates the number of occurences of b

We use vec(a*(b ⊗ b ⊗ ... ⊗ b)) = (b' ⊗ b' ⊗ ... ⊗ b' ⊗ I)vec(a)

"""
function a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
    ma, na = size(a)
    mb, nb = size(b)
    mc, nc = size(c)
    mborder = mb^order
    nborder = nb^order
    na == mborder || throw(DimensionMismatch("The number of columns of a, $na, doesn't match the number of rows of b, $mb, times order = $order"))
    mc == ma || throw(DimensionMismatch("The number of rows of c, $mc, doesn't match the number of rows of a, $ma"))
    nc == nborder || throw(DimensionMismatch("The number of columns of c, $nc, doesn't match the number of columns of b, $nb, times order = $order"))

    # copy input
    copyto!(work1,a)
    n1 = ma*na
    n2 = Int(n1*nb/mb)
    v1 = view(work1,1:n1)
    v2 = view(work2,1:n2)
    for q=0:order-1
        kron_mul_elem_t!(v2,b,v1,mb^(order-q-1),nb^q*ma)
        if q < order - 1
            # update dimension
            n2 = Int(n2*nb/mb)
            # swap and resize work vectors
            v1parent = v1.parent
            v1 = v2
            v2 = view(v1parent,1:n2)
        end
    end
    copy!(c,v2)
end

function a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, order::Int64)
    ma, na = size(a)
    mb, nb = size(b)
    mc, nc = size(c)
    mborder = mb^order
    nborder = nb^order
    na == mborder || throw(DimensionMismatch("The number of columns of a, $na, doesn't match the number of rows of b, $mb, times order = $order"))
    mc == ma || throw(DimensionMismatch("The number of rows of c, $mc, doesn't match the number of rows of a, $ma"))
    nc == nborder || throw(DimensionMismatch("The number of columns of c, $nc, doesn't match the number of columns of b, $nb, times order = $order"))
    mb == nb || throw(DimensionMismatch("B must be a square matrix"))

    avec = vec(a)
    cvec = vec(c)
    for q=0:order-1
#        vavec = view(avec,1:ma*mb^(order-q)*nb^q)
#        vcvec = view(cvec,1:ma*mb^(order-q-1)*nb^(q+1))
        kron_mul_elem_t!(cvec,b,avec,mb^(order-q-1),nb^q*ma)
        if q < order - 1
            copy!(avec,cvec)
        end
    end
end

"""
    a_mul_b_kron_c!(d::AbstractVecOrMat, a::AbstractVecOrMat, b::AbstractMatrix, c::AbstractMatrix, order::Int64)

Performs a*B*(c ⊗ c ⊗ ... ⊗ c). The solution is returned in matrix or vector d. order indicates the number of occurences of c. c must be a square matrix

We use vec(a*b*(c ⊗ c ⊗ ... ⊗ c)) = (c' ⊗ c' ⊗ ... ⊗ c' ⊗ a)vec(b)

"""
function a_mul_b_kron_c!(d::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, order::Int64,work1::AbstractVector,work2::AbstractVector)
    ma, na = size(a)
    mb, nb = size(b)
    mc, nc = size(c)
    md, nd = size(d)
    na == mb || throw(DimensionMismatch("The number of columns of a, $(size(a,2)), doesn't match the number of rows of b, $(size(b,1))"))
    nb == mc^order || throw(DimensionMismatch("The number of columns of b, $(size(b,2)), doesn't match the number of rows of c, $(size(c,1)), times order, $order"))
    (ma == md && nc^order == nd) || throw(DimensionMismatch("Dimension mismatch for D: $(size(d)) while ($ma, $(nc^order)) was expected"))
    mul!(reshape(view(work1,1:ma*nb),ma,nb),a,b)
    for q=0:order-1
        kron_mul_elem_t!(work2,c,work1,mc^(order-q-1),nc^q*ma)
        if q < order - 1
            copy!(work1,work2)
        end
    end
    copyto!(d,work2)
end

"""
    function kron_at_kron_b_mul_c!(d::AbstractVector, offset_d::Int64, a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, offset_c::Int64, work1::AbstractVector, work2::AbstractVector, offset_w::Int64)
computes d = (a^T ⊗ a^T ⊗ ... ⊗ a^T ⊗ b)c using data from c at offset_c and work vectors work1 and work2 at offset_w
"""
function kron_at_kron_b_mul_c!(d::AbstractVector, offset_d::Int64, a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, offset_c::Int64, work1::AbstractVector, work2::AbstractVector, offset_w::Int64)
    mb,nb = size(b)
    if order == 0
        A_mul_B!(d,offset_d,b,1,mb,nb,c,offset_c,1)
    else
        ma, na = size(a)
        #    length(work) == naorder*mb  || throw(DimensionMismatch("The dimension of vector , $(length(c)) doesn't correspond to order, ($order)  and the dimension of the matrices a, $(size(a)), and b, $(size(b))"))
        p = ma^order
        kron_mul_elem!(work1, 1, b, c, offset_c, p, 1)
        p = Int(p/ma)
        q = mb
        for i = 1:order
            kron_mul_elem_t!(work2, offset_w, a, work1, 1, p, q)
            if i < order
                copyto!(work1, 1, work2, offset_w, p*na*q)
                p = Int(p/ma)
                q *= na
            end
        end
        copyto!(d, offset_d, work2, offset_w, p*na*q)
    end
end

"""
    function kron_at_kron_b_mul_c!(d::AbstractVector, a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, work1::AbstractVector, work2::AbstractVector)
computes d = (a^T ⊗ a^T ⊗ ... ⊗ a^T ⊗ b)c using work vectors work1 and work2
"""
function kron_at_kron_b_mul_c!(d::AbstractVector, a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, work1::AbstractVector, work2::AbstractVector)
    mb,nb = size(b)
    if order == 0
        A_mul_B!(d,1,b,1,mb,nb,c,1,1)
    else
        ma, na = size(a)
        #    length(work) == naorder*mb  || throw(DimensionMismatch("The dimension of vector , $(length(c)) doesn't correspond to order, ($order)  and the dimension of the matrices a, $(size(a)), and b, $(size(b))"))
        p = ma^order
        kron_mul_elem!(work1, b, c, p, 1)
        p = Int(p/ma)
        q = mb
        for i = 1:order
            kron_mul_elem_t!(work2, a, work1, p, q)
            if i < order
                copyto!(work1, 1, work2, 1, p*na*q)
                p = Int(p/ma)
                q *= na
            end
        end
        copyto!(d,1,work2,1,p*na*q)
    end
end

"""
    function kron_at_kron_b_mul_c!(a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, offset_c::Int64, work::AbstractVector)
updates c at offset_c with (a^T ⊗ a^T ⊗ ... ⊗ a^T ⊗ b)c using c and work as work vectors
"""
function kron_at_kron_b_mul_c!(a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, offset_c::Int64, work::AbstractVector)
    kron_at_kron_b_mul_c!(work, 1, a, order, b, c, offset_c, work, c, offset_c)
    copy!(c, work)
end

"""
    function kron_at_kron_b_mul_c!(a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, work::AbstractVector)
updates c  with (a^T ⊗ a^T ⊗ ... ⊗ a^T ⊗ b)c using c and work as work vectors
"""
function kron_at_kron_b_mul_c!(a::AbstractMatrix, order::Int64, b::AbstractMatrix, c::AbstractVector, work::AbstractVector)
    kron_at_kron_b_mul_c!(work, a, order, b, c, work, c)
    copy!(c, work)
end

function kron_at_mul_b!(c::AbstractVector, a::AbstractMatrix, order::Int64, b::AbstractVector, q::Int64, work1::AbstractVector, work2::AbstractVector)
    ma,na = size(a)
#    length(work) == naorder*mb  || throw(DimensionMismatch("The dimension of vector , $(length(c)) doesn't correspond to order, ($order)  and the dimension of the matrices a, $(size(a)), and b, $(size(b))"))
    p = ma^(order-1)
    s = p*q*na
    copyto!(work1,b)
    for i = 1:order
        kron_mul_elem_t!(work2, a, work1, p, q)
        if i < order
            copyto!(work1, 1, work2, 1, s)
            p = Int(p/ma)
            q *= na
            s = p*q*na
        end
    end
    copyto!(c,1,work2,1,s)
end

function kron_a_mul_b!(c::AbstractVector, a::AbstractMatrix, order::Int64, b::AbstractVector, q::Int64, work1::AbstractVector, work2::AbstractVector)
    ma,na = size(a)
#    length(work) == naorder*mb  || throw(DimensionMismatch("The dimension of vector , $(length(c)) doesn't correspond to order, ($order)  and the dimension of the matrices a, $(size(a)), and b, $(size(b))"))
    p = na^(order-1)
    s = p*q*ma
    copyto!(work1,b)
    for i = 1:order
        kron_mul_elem!(work2,a,work1,p,q)
        if i < order
            copyto!(work1, 1, work2, 1, s)
            p = Int(p/na)
            q *= ma
            s = p*q*ma
        end
    end
    copyto!(c,1,work2,1,s)
end

function at_mul_b_kron_c!(d::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
    ma, na = size(a)
    mb, nb = size(b)
    mc, nc = size(c)
    if mc <= nc
        At_mul_B!(work1, 1, a, 1, na, ma, b, 1, nb)
        kron_at_mul_b!(vec(d), c, order, work1, na, vec(d), work2)
    else
        kron_at_mul_b!(work1, c, order, b, na, work1, work2)
        At_mul_B!(vec(d), 1, a, 1, na, ma, work1, 1, nc^order)
    end
end

function a_mul_b_kron_ct!(d::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)
    ma, na = size(a)
    mb, nb = size(b)
    mc, nc = size(c)
    if mc <= nc
        A_mul_B!(work1, 1, a, 1, ma, na, b, 1, nb)
        kron_a_mul_b!(vec(d), c, order, work1, ma, work1, work2)
    else
        kron_a_mul_b!(work1, b, order, c, mb, work1, work2)
        At_mul_B!(vec(d), 1, a, 1, na, ma, work1, 1, nc^order)
    end
end

function a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::AbstractVector,work::AbstractVector)
    ma, na = size(a)
    order = length(b)
    mc, nc = size(c)
    mborder = 1
    nborder = 1
    for i = 1:order
        mb, nb = size(b[i])
        mborder *= mb
        nborder *= nb
    end
    na == mborder || throw(DimensionMismatch("The number of columns of a, $na, doesn't match the number of rows of matrices in b, $mborder"))
    mc == ma || throw(DimensionMismatch("The number of rows of c, $mc, doesn't match the number of rows of a, $ma"))
    nc == nborder || throw(DimensionMismatch("The number of columns of c, $nc, doesn't match the number of columns of matrices in b, $nborder"))
    mborder <= nborder || throw(DimensionMismatch("the product of the number of rows of the b matrices needs to be smaller or equal to the product of the number of columns. Otherwise, you need to call a_mul_kron_b!(c::AbstractMatrix, a::AbstractMatrix, b::Vector{AbstractMatrix},work::AbstractVector)"))

    copyto!(work,a)
    mb, nb = size(b[1])
    vwork = view(work,1:ma*na)
    cvec = vec(c)
    mc = ma*na
    p = Int(mborder/size(b[order],1))
    q = ma
    for i = order:-1:1
        mb, nb = size(b[i])
        mc = Int(mc*nb/mb)
        vcvec = view(cvec,1:mc)
        kron_mul_elem_t!(cvec,b[i],work,p,q)
        if i > 1
            p = Int(p/mb)
            q = q*nb
            vwork = view(work,1:mc)
            copy!(vwork,vcvec)
        end
    end
end

"""
    a_mul_b_kron_c_d!(d::AbstractVecOrMat, a::AbstractVecOrMat, b::AbstractMatrix, c::AbstractMatrix, order::Int64)

Performs e = a*b*(c ⊗ d ⊗ ... ⊗ d). The solution is returned in matrix or vector e order indicates the number of occurences of d. c and d must be square matrices

We use vec(a*b*(c ⊗ d ⊗ ... ⊗ d)) = (c' ⊗ d' ⊗ ... ⊗ d' ⊗ a)vec(b)

"""
function a_mul_b_kron_c_d!(e::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::AbstractMatrix, d::AbstractMatrix, order::Int64, work1::AbstractVector, work2::AbstractVector)

     ma, na = size(a)
     mb, nb = size(b)
     mc, nc = size(c)
     md, nd = size(d)
     me, ne = size(e)
     na == mb || throw(DimensionMismatch("The number of columns of a, $(size(a,2)), doesn't match the number of rows of b, $(size(b,1))"))
     nb == mc*md^(order-1) || throw(DimensionMismatch("The number of columns of b, $(size(b,2)), doesn't match the number of rows of c, $(size(c,1)), and d, $(size(d,1)) for order, $order"))
     (ma == me && nc*nd^(order-1) == ne) || throw(DimensionMismatch("Dimension mismatch for e: $(size(e)) while ($ma, $(nc*nd^(order-1))) was expected"))
     A_mul_B!(work1, 1, a, 1, ma, mb, vec(b), 1, nb)
     p = mc*md^(order - 2)
     q = ma
     for i = 0:order - 2
        kron_mul_elem_t!(work2, 1, d, work1, 1, p, q)
        copy!(work1,work2)
        p = Int(p/md)
        q *= nd
    end
    kron_mul_elem_t!(work2, c, work1, 1,q)
    copyto!(e, 1, work2, 1, ma*nc*nd^(order-1))
end

convert(::Type{Array{Float64, 2}}, x::Base.ReshapedArray{Float64,2,SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true},Tuple{}}) = unsafe_wrap(Array,pointer(x.parent.parent,x.parent.indexes[1][1]),x.dims)

end


text:
type Controller
    ctx::Context
    sock::Socket
    port::Int
    key::String
end

function Controller(port)
    ctx = Context()
    sock = Socket(ctx, REQ)
    connect(sock, "tcp://localhost:$port")
    key = create_key(sock)
    return Controller(ctx, sock, port, key)
end

Base.show(io::IO, ctr::Controller) =
    print(io, "Controller($(ctr.sock.data),$(ctr.port))")

@doc """Send control message""" ->
function control(sock::Socket, msg::String)
    send(sock, "c:" * msg)
    resp = bytestring(recv(sock))
    return resp
end

@doc """Send control message""" ->
function control(ctr::Controller, msg::String)
    return control(ctr.sock, msg)
end

@doc """Create new buffer on Collector, return its key""" ->
function create_key(sock::Socket)
    key = randstring(5)
    control(sock, "createkey $key")
    return key
end

@doc """Get all data and remove buffers on Collector""" ->
function finalize(ctr::Controller)
    send(ctr.sock, "c:finalize $(ctr.key)")
    # read all data from socket
    stop = false
    timetable = String[]
    while !stop
        line = bytestring(recv(ctr.sock))
        if line == "--END--"
            stop = true
        else
            push!(timetable, line)
        end
        send(ctr.sock, "ok")
    end
    return timetable
end


text:
using Avro
using Test

include("schema.jl")
include("io.jl")
include("generic.jl")
include("data_file.jl")


text:
function doji(ohlc::TimeArray; width = .01, o = :Open, h = :High, l = :Low, c = :Close)
  res = @. abs((ohlc[o] - ohlc[c])) / (ohlc[h] - ohlc[l]) < width
  TimeArray(timestamp(res), values(res), [:doji], meta(ohlc))
end

function hammer(x)
  #code here
end

function inverted_hammer(x)
  #code here
end

function hanging_man(x)
  #code here
end

function mirabosu(x)
  #code here
end

function shooting_star(x)
  #code here
end

function spinning_top(x)
  #code here
end

function three_white_soldiers(x)
  #code here
end

function three_black_crows(x)
  #code here
end

function morning_star(x)
  #code here
end


text:
###### GOAL: #######
### By exchanging bits of 2 bitvectors at the same position,
### minimize the number of breakpoints in both strands, and
### ruturn that number.
####################
# Example 1
# x = [0 1 1 1]
# y = [1 0 0 0]
# We can exchange position 1:
# x = [1 1 1 1]
# y = [0 0 0 0]
# then breakpoint = 0
#
# Example 2
# x = [0 1 0 0 1]
# y = [1 0 1 0 0]
# we can exchange position 2:
# x = [0 0 0 0 1]
# y = [1 1 1 0 0]
# then breakpoint = 2 (1 at end of x and another at middle of y)

using Random
using Test
using BenchmarkTools

"""
Starting at position `n`, search through to the end of `seq1` and `seq2` to see if
position `seq1[n - 1]` and `seq2[n - 1]` should be swapped. Swapping index are
stored in `flip_idx`.

Returns the total number of breakpoints in `seq1` and `seq2` after swapping.
"""
function binary_flip!(n, seq1, seq2, flip_idx,
    intermediates::Vector{Int} = [-1 for i in 1:(length(seq1) + 1)]
    )
    # quick lookup
    if intermediates[n] != -1
        return intermediates[n]
    end

    if n > length(seq1)
        return 0
    elseif n == 1
        # only flip bits in previous position
        return binary_flip!(n + 1, seq1, seq2, intermediates)
    elseif seq1[n] == seq2[n]
        # don't flip and only calculate error
        return binary_flip!(n + 1, seq1, seq2, flip_idx, intermediates) + !(seq1[n] == seq1[n - 1]) + !(seq2[n] == seq2[n - 1])
    else
        # calculate error of flip/noflip by recursion
        yesflip = binary_flip!(n + 1, seq1, seq2, flip_idx, intermediates) + !(seq1[n] == seq2[n - 1]) + !(seq2[n] == seq1[n - 1])
        noflip  = binary_flip!(n + 1, seq1, seq2, flip_idx, intermediates) + !(seq1[n] == seq1[n - 1]) + !(seq2[n] == seq2[n - 1])
        # println("n = $n, yesflip = $yesflip, noflip = $noflip, min = $(min(yesflip, noflip))")
        # store intermediate results for later retrival
        intermediates[n] = min(yesflip, noflip)
        if yesflip < noflip
            # record flipping location, flip 2 sequence at previous location,
            flip_idx[n - 1] = true
            seq1[n - 1], seq2[n - 1] = seq2[n - 1], seq1[n - 1]
            # println("flipped $(n - 1)!")
        end
        return min(yesflip, noflip)
    end
end
binary_flip!(seq1, seq2, flip) = binary_flip!(2, seq1, seq2, flip) #start at position 2

x = [1 0]
y = [0 1]
flip = falses(2)
@test binary_flip!(x, y, flip) == 0
@test flip == [true; false]
@test x == [0 0]
@test y == [1 1]

x = [1 0 1]
y = [0 1 0]
flip = falses(3)
@test binary_flip!(x, y, flip) == 0
@test flip == [false; true; false]
@test x == [1 1 1]
@test y == [0 0 0]

x = [0 1 1 1]
y = [1 0 0 0]
flip = falses(4)
@test binary_flip!(x, y, flip) == 0
@test flip == [true; false; false; false]
@test x == [1 1 1 1]
@test y == [0 0 0 0]

x = [1 0 1 0]
y = [0 1 0 0]
flip = falses(4)
@test binary_flip!(x, y, flip) == 1
@test flip == [false; true; false; false]
@test x == [1 1 1 0]
@test y == [0 0 0 0]

x = [0 1 0 0 1]
y = [1 0 1 0 0]
flip = falses(5)
@test binary_flip!(x, y, flip) == 2
@test flip == [false; true; false; false; false]
@test x == [0 0 0 0 1]
@test y == [1 1 1 0 0]

x = [1 0 1 0]
y = [0 1 0 1]
flip = falses(4)
@test binary_flip!(x, y, flip) == 0
@test flip == [true; false; true; false]
@test x == [0 0 0 0]
@test y == [1 1 1 1]

x = [0 1 1 1 0 0]
y = [1 1 1 0 0 0]
flip = falses(6)
@test binary_flip!(x, y, flip) == 3
@test flip == [false; false; false; false; false; false]
@test x == [0 1 1 1 0 0]
@test y == [1 1 1 0 0 0]

x = [0 1 0 1 1 0]
y = [1 0 1 0 0 0]
flip = falses(6)
@test binary_flip!(x, y, flip) == 1
@test flip == [true; false; true; false; false; false]
@test x == [1 1 1 1 1 0]
@test y == [0 0 0 0 0 0]

x = [0 1 1 1 1 0]
y = [1 0 1 0 0 0]
flip = falses(6)
@test binary_flip!(x, y, flip) == 3
@test flip == [true; false; false; false; false; false]
@test x == [1 1 1 1 1 0]
@test y == [0 0 1 0 0 0]

x = [0 1 1 1 1 0]
y = [1 0 1 0 0 1]
flip = falses(6)
@test binary_flip!(x, y, flip) == 2
@test flip == [true; false; false; true; true; false]
@test x == [1 1 1 0 0 0]
@test y == [0 0 1 1 1 1]

x = [0 1 1 1 1 0]
y = [1 0 1 1 0 1]
flip = falses(6)
@test binary_flip!(x, y, flip) == 2
@test flip == [true; false; false; false; true; false]
@test x == [1 1 1 1 0 0]
@test y == [0 0 1 1 1 1]

# nonmemoized (bitrand(100) would fail)
x = bitrand(30)
y = bitrand(30)
flip = falses(30)
@btime binary_flip!($x, $y, $flip) # 298.914 μs (1 allocation: 336 bytes)

# memoized
x = bitrand(30)
y = bitrand(30)
@btime binary_flip!($x, $y, $flip) # 737.567 ns (1 allocation: 336 bytes)

x = bitrand(10000)
y = bitrand(10000)
@time binary_flip!(x, y)

Random.seed!(2020)
x = bitrand(10)
y = bitrand(10)
[x y]
flip = falses(10);
@time binary_flip!(x, y, flip)
[x y]

Random.seed!(2020)
x = bitrand(20)
y = bitrand(20)
[x y]
flip = falses(20);
@time binary_flip!(x, y, flip)
[x y]

#more examples
x = bitrand(10)
y = bitrand(10)
[x y]
flip = falses(10);
@time binary_flip!(x, y, flip)
[x y]

x = bitrand(10000)
y = bitrand(10000)
[x y]
flip = falses(10000);
@time binary_flip!(x, y, flip)
[x y]

function a()
    println("reached a!")
    return 1
end
function b()
    println("reached b!")
    return 2
end
c = a() + b()
# right hand side of equations are evaluated from left to right
# julia> c = a() + b()
# reached a!
# reached b!
# 3

function has_intersect!(c::BitVector, a::BitVector, b::BitVector)
    c .= a .& b
    return any(c)
end

function has_intersect(a::BitVector, b::BitVector)
    @inbounds for i in eachindex(a)
        if a[i] && b[i]
            return true
        end
    end
    return false
end

a = falses(100_000)
b = falses(100_000)
c = falses(100_000)

Random.seed!(2020)
a[1:1000] .= b[1:1000] .= true
shuffle!(a)
shuffle!(b)

@btime has_intersect($a, $b)
@btime has_intersect!($c, $a, $b)

function set_flip!(
    n::Int,
    vector_set1::Vector{BitVector},
    vector_set2::Vector{BitVector},
    flip_idx::BitVector,
    intermediates::Vector{Int} = [-1 for i in 1:(length(vector_set1) + 1)],
    storage::BitVector = falses(length(vector_set1[1]))
    )
    # quick lookup
    if intermediates[n] != -1
        return intermediates[n]
    end

    if n > length(vector_set1)
        return 0
    elseif n == 1
        # only flip bits in previous position
        return set_flip!(n + 1, vector_set1, vector_set2, flip_idx, intermediates)
    elseif vector_set1[n] == vector_set2[n]
        # don't flip and only calculate error
        return set_flip!(n + 1, vector_set1, vector_set2, flip_idx, intermediates) +
                !has_intersect!(storage, vector_set1[n], vector_set1[n - 1]) +
                !has_intersect!(storage, vector_set2[n], vector_set2[n - 1])
    else
        # calculate error of flip/noflip by recursion
        yesflip = set_flip!(n + 1, vector_set1, vector_set2, flip_idx, intermediates) +
                    !has_intersect!(storage, vector_set1[n], vector_set2[n - 1]) +
                    !has_intersect!(storage, vector_set2[n], vector_set1[n - 1])
        noflip  = set_flip!(n + 1, vector_set1, vector_set2, flip_idx, intermediates) +
                    !has_intersect!(storage, vector_set1[n], vector_set1[n - 1]) +
                    !has_intersect!(storage, vector_set2[n], vector_set2[n - 1])

        # store intermediate results for later retrival
        intermediates[n] = min(yesflip, noflip)
        if yesflip < noflip
            # record flipping location, flip 2 sequence at previous location,
            flip_idx[n - 1] = true
            vector_set1[n - 1], vector_set2[n - 1] = vector_set2[n - 1], vector_set1[n - 1]
        end
        return min(yesflip, noflip)
    end
end
set_flip!(vector_set1, vector_set2, flip) = set_flip!(2, vector_set1, vector_set2, flip) #start at position 2

function has_intersect!(c::BitVector, a::BitVector, b::BitVector)
    c .= a .& b
    return any(c)
end

using Random
using Test
using BenchmarkTools

Random.seed!(2020)

x = [bitrand(10), bitrand(10), bitrand(10), bitrand(10)]
y = [bitrand(10), bitrand(10), bitrand(10), bitrand(10)]
flip = falses(4);
@test set_flip!(x, y, flip) == 2
@test flip == [false; true; false; false]

x = [bitrand(5), bitrand(5), bitrand(5), bitrand(5), bitrand(5), bitrand(5), bitrand(5)]
y = [bitrand(5), bitrand(5), bitrand(5), bitrand(5), bitrand(5), bitrand(5), bitrand(5)]
flip = falses(7);
@test set_flip!(x, y, flip) == 1
@test flip == [false; false; false; false; true; false; false]

using Random
using Test
using BenchmarkTools

Random.seed!(1)

a = bitrand(5)
b = bitrand(5)

struct bitvec_pair

end

for pair in (a, b)
    println(pair)
end

max_stars_repo_name: Hua-Zhou/MendelImpute


text:
# Currently supported distributions
const DISTS = [
    Beta,
    BetaFourParameters,
    Exponential,
    Gamma,
    LogitNormal,
    LogNormal,
    NegativeBinomial,
    Normal,
    Poisson,
    TDist,
    TDistLocationScale,
    Weibull
]

export Beta,
    BetaFourParameters,
    Exponential,
    Gamma,
    LogitNormal,
    LogNormal,
    NegativeBinomial,
    Normal,
    Poisson,
    TDist,
    TDistLocationScale,
    Weibull

"""
    score!(score_til::Matrix{T}, y::T, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T

Fill `score_til` with the score of distribution `D` with parameters `param[:, t]` considering
the observation `y`.
"""
function score!(score_til::Matrix{T}, y::T, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
    return error("score! not implemented for $D distribution")
end

function score!(score_til::Matrix{T}, y::Int, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
    return error("score! not implemented for $D distribution")
end

"""
    fisher_information!(aux::AuxiliaryLinAlg{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T

Fill `aux` with the fisher information of distribution `D` with parameters `param[:, t]`.
"""
function fisher_information!(aux::AuxiliaryLinAlg{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
    return error("fisher_information! not implemented for $D distribution")
end

"""
    log_likelihood(D::Type{<:Distribution}, y::Vector{T}, param::Matrix{T}, n::Int) where T

Evaluate the log-likelihood of the distribution `D` considering the time-varying parameters
`param` and the observations `y`.
"""
function log_likelihood(D::Type{<:Distribution}, y::Vector{T}, param::Matrix{T}, n::Int) where T
    return error("log_likelihood not implemented for $D distribution")
end

"""
    link!(param_tilde::Matrix{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T

Fill `param_tilde` after the unlinking procedure of `param`.
"""
function link!(param_tilde::Matrix{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
    return error("link! not implemented for $D distribution")
end

"""
    unlink!(param::Matrix{T}, D::Type{<:Distribution}, param_tilde::Matrix{T}, t::Int) where T

Fill `param` after the unlinking procedure of `param_tilde`.
"""
function unlink!(param::Matrix{T}, D::Type{<:Distribution}, param_tilde::Matrix{T}, t::Int) where T
    return error("unlink! not implemented for $D distribution")
end

"""
    jacobian_link!(aux::AuxiliaryLinAlg{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T

Write the jacobian of the link map in `aux`.
"""
function jacobian_link!(aux::AuxiliaryLinAlg{T}, D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
    return error("jacobian_link! not implemented for $D distribution")
end

"""
    update_dist(D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T

Create a new distribution from Distributions.jl based on the parametrization used in
ScoreDrivenModels.jl.
"""
function update_dist(D::Type{<:Distribution}, param::Matrix{T}, t::Int) where T
    return error("update_dist not implemented for $D distribution")
end

"""
    params_sdm(d::Distribution)

Recover the parametrization used in ScoreDrivenModels.jl based on a Distribution from
Distributions.jl.
"""
function params_sdm(d::Distribution)
    return error("params_sdm not implemented for $d distribution")
end

"""
    num_params(D::Type{<:Distribution})

Number of parameters of a given distribution.
"""
function num_params(D::Type{<:Distribution})
    return error("num_params not implemented for $D distribution")
end


text:
module Lupin

export partition
include("partition.jl")

export vecjoin
include("vecjoin.jl")

include("mergewith.jl")

include("dict.jl")

include("compose.jl")

include("enums.jl")

include("reverse.jl")

export Argument, apply
include("arguments.jl")

end # module Lupin