text:
"""
    pmkstruct(P; CF1, grade=l, fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> KRInfo, iz, ip

Determine the Kronecker-structure and infinite pole-zero structure information of the polynomial matrix `P(λ)` and
return an `KRInfo` object and the multiplicities of infinite zeros and poles.
The computation of the Kronecker-structure employs strong linearizations of `P(λ)` in either
the first companion form, if `CF1 = true`, or the second companion form, if `CF1 = false`.
The effective grade `l` to be used for linearization can be specified via the keyword argument `grade` as
`grade = l`, where `l` must be chosen equal to or greater than the degree of `P(λ)`.
The default value used for `l` is `l = deg(P(λ))`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The information on the Kronecker-structure consists of the right Kronecker indices `rki`,
left Kronecker indices `lki`, infinite elementary divisors `id`, the
number of finite eigenvalues `nf` and normal rank `nrank` and can be obtained from `KRInfo` as
`KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`KRInfo`](@ref).

The determination of the Kronecker-structure information is performed by building a companion
form linearization `M-λN` of `P(λ)` and reducing the pencil `M-λN` to an
appropriate Kronecker-like form (KLF) which exhibits the number of finite eigenvalues,
the multiplicities of the infinite eigenvalues,
the left and right Kronecker indices and the normal rank.
The Kronecker-structure information and pole-zero stucture information on `P(λ)` are recovered from
the Kronecker-structure information on `M-λN` using the results of [1] and [2].

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

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

The multiplicities of infinite eigenvalues are provided in the integer vector `id`,
where each `i`-th element `id[i]` is the order of an infinite elementary divisor
(i.e., the multiplicity of an infinite eigenvalue).

The multiplicities of the infinite zeros of `P(λ)` are returned in `iz` and represent
the positive differences between the multiplicities of
the infinite eigenvalues of `P(λ)` and the effective grade `l` of `P(λ)` [2].

The multiplicities of the infinite poles of `P(λ)` are returned in `ip` and represent
the absolute values of the negative differences between the multiplicities of
the infinite eigenvalues of `P(λ)` and the effective grade `l` of `P(λ)` [2].

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

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

[1] F. De Terán, F. M. Dopico, D. S. Mackey, Spectral equivalence of polynomial matrices and
the Index Sum Theorem, Linear Algebra and Its Applications, vol. 459, pp. 264-333, 2014.

[2] A. Varga, On computing the Kronecker structure of polynomial matrices using Julia, June 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function pmkstruct(P::AbstractArray{T,3}; CF1::Bool = size(P,1) <= size(P,2) ? false : true,
                   grade::Int = pmdeg(P), fast::Bool = false, atol::Real = zero(real(T)),
                   rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   if CF1
      M, N = pm2lpCF1(P, grade = grade)
      kinfo = pkstruct(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      kinfo.rki == [] || (grade > 1 && (kinfo.rki .-= grade-1))
      kinfo.nrank > 0 && grade > 1 && (kinfo.nrank -= size(P,2)*(grade-1))
   else
      M, N = pm2lpCF2(P, grade = grade)
      kinfo = pkstruct(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      kinfo.lki == [] || (grade > 1 && (kinfo.lki .-= grade-1))
      kinfo.nrank > 0 && grade > 1 && (kinfo.nrank -= size(P,1)*(grade-1))
   end
   # determine the infinite structural indices
   grade < 0 ? is = Int[] : (is = [zeros(Int,kinfo.nrank-length(kinfo.id)); kinfo.id] .- grade)
   iz = is[is.>0]
   ip = sort(abs.(is[is.<0]))

   return kinfo, iz, ip
end
pmkstruct(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
          pmkstruct(poly2pm(P); kwargs...)
"""
    pmeigvals(P; CF1, grade = l, fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, KRInfo)

Return the finite and infinite eigenvalues of the polynomial matrix `P(λ)` in `val`
and information on the Kronecker-structure of `P(λ)` in the `KRInfo` object.
The computation of eigenvalues and Kronecker-structure employs strong linearizations of `P(λ)` in either
the first companion form, if `CF1 = true`, or the second companion form, if `CF1 = false`.
The effective grade `l` to be used for linearization can be specified via the keyword argument `grade` as
`grade = l`, where `l` must be chosen equal to or greater than the degree of `P(λ)`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The information on the Kronecker-structure consists of the right Kronecker indices `rki`,
left Kronecker indices `lki`, infinite elementary divisors `id`, the
number of finite eigenvalues `nf` and normal rank `nrank` and can be obtained from `KRInfo` as
`KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively. For more details, see  [`KRInfo`](@ref).

The computation of the eigenvalues is performed by building the companion form based linearization `M-λN`
of the polynomial matrix `P(λ)` and then reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF)
which exhibits the number of finite eigenvalues, the multiplicities of the infinite eigenvalues,
the left and right Kronecker indices and the normal rank.
The left and right Kronecker indices of `P(λ)` are returned in `KRInfo.rki` and `KRInfo.rki`, respectively, and
their values are recovered from the left and right Kronecker indices of `M-λN` using the results of [1].
The multiplicities of the infinite eigenvalues of `P(λ)` of effective grade `l` is returned in `KRInfo.id`.
The number of finite eigenvalues in `val` is equal to the number of finite eigenvalues of `M-λN` (returned in `KRInfo.nf`),
while the number of infinite eigenvalues in `val` is the sum of multiplicites returned in `KRInfo.id`.

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

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

[1] F. De Terán, F. M. Dopico, D. S. Mackey, Spectral equivalence of polynomial matrices and
the Index Sum Theorem, Linear Algebra and Its Applications, vol. 459, pp. 264-333, 2014.
"""
function pmeigvals(P::AbstractArray{T,3}; CF1::Bool = size(P,1) <= size(P,2) ? false : true,
                   grade::Int = pmdeg(P), fast::Bool = false, atol::Real = zero(real(T)),
                   rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   if CF1
      M, N = pm2lpCF1(P, grade = grade)
      val, kinfo = peigvals(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      kinfo.rki == [] || (grade > 1 && (kinfo.rki .-= grade-1))
      kinfo.nrank > 0 && grade > 1 && (kinfo.nrank -= size(P,2)*(grade-1))
   else
      M, N = pm2lpCF2(P, grade = grade)
      val, kinfo = peigvals(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      kinfo.lki == [] || (grade > 1 && (kinfo.lki .-= grade-1))
      kinfo.nrank > 0 && grade > 1 && (kinfo.nrank -= size(P,1)*(grade-1))
   end
   return val, kinfo
end
pmeigvals(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
          pmeigvals(poly2pm(P); kwargs...)
"""
    pmzeros(P; CF1, fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the finite and infinite zeros of the polynomial matrix `P(λ)` in `val`,
the multiplicities of infinite zeros in `iz` and
information on the  Kronecker-structure of `P(λ)` in the `KRInfo` object.
The computation of zeros and Kronecker-structure employs strong linearizations of `P(λ)` in either
the first companion form, if `CF1 = true`, or the second companion form, if `CF1 = false`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The information on the Kronecker-structure consists of the right Kronecker indices `rki`,
left Kronecker indices `lki`, infinite elementary divisors `id`, the
number of finite eigenvalues `nf`, the normal rank `nrank`, and can be obtained from `KRInfo` as
`KRInfo.rki`, `KRInfo.lki`, `KRInfo.id`, `KRInfo.nf` and `KRInfo.nrank`, respectively.
For more details, see  [`KRInfo`](@ref).

The computation of the zeros is performed by building the companion form based linearization `M-λN`
of the polynomial matrix `P(λ)` and then reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF)
which exhibits information on its Kronecker structure.
The left and right Kronecker indices of `P(λ)` are returned in `KRInfo.rki` and `KRInfo.rki`, respectively, and
their values are recovered from the left and right Kronecker indices of `M-λN` using the results of [1].
The multiplicities of the infinite zeros of `P(λ)`, returned in `iz`, are the positive differences between the multiplicities of
the infinite eigenvalues of `M-λN` and the degree of `P(λ)` [2].
The number of finite zeros in `val` is equal to the number of finite eigenvalues of `M-λN` (returned in `KRInfo.nf`),
while the number of infinite zeros in `val` is the sum of multiplicites in `iz`.

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

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

[1] F. De Terán, F. M. Dopico, D. S. Mackey, Spectral equivalence of polynomial matrices and
the Index Sum Theorem, Linear Algebra and Its Applications, vol. 459, pp. 264-333, 2014.

[2] A. Varga, On computing the Kronecker structure of polynomial matrices using Julia, June 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function pmzeros(P::AbstractArray{T,3}; CF1::Bool = size(P,1) <= size(P,2) ? false : true,
                 fast = false, atol::Real = zero(real(T)),
                 rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   dm1 = pmdeg(P)-1
   if CF1
      M, N = pm2lpCF1(P)
      val, iz, kinfo = pzeros(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      kinfo.rki == [] || dm1 <= 0 || (kinfo.rki .-= dm1)
      kinfo.nrank > 0 && dm1 > 0 && (kinfo.nrank -= size(P,2)*dm1)
   else
      M, N = pm2lpCF2(P)
      val, iz, kinfo = pzeros(M, N, fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
      kinfo.lki == [] || dm1 <= 0 || (kinfo.lki .-= dm1)
      kinfo.nrank > 0 && dm1 > 0 && (kinfo.nrank -= size(P,1)*dm1)
   end
   if iz != [] && dm1 > 0
      iz = iz[iz .> dm1] .- (dm1)
      val = [val[1:kinfo.nf]; Inf*ones(real(T),sum(iz))]
   end
   return val, iz, kinfo
end
pmzeros(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       pmzeros(poly2pm(P); kwargs...)
"""
    pmpoles(P; CF1, fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, ip, id)

Return the finite and infinite poles of the polynomial matrix `P(λ)` in `val`,
the multiplicities of infinite poles in `ip` and the infinite elementary divisors of
`Q(λ) = [P(λ) I; I 0]` in `id`.
The computation of pole-structure employs strong linearizations of `Q(λ)` in either
the first companion form, if `CF1 = true`, or the second companion form, if `CF1 = false`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The information on the finite and infinite poles of the polynomial matrix `P(λ)` is obtained from the
Kronecker-structure information of the underlying linearization `Q(λ)` using the results of [1] and [2].
The computation of the poles is performed by building the companion form based linearization `M-λN`
of the extended polynomial matrix `Q(λ) = [P(λ) I; I 0]` (see [1]) and then reducing the pencil `M-λN` to an appropriate
Kronecker-like form (KLF) which exhibits information on its Kronecker structure.
Since `Q(λ)` is regular, `M-λN` has no left or right Kronecker indices.
The multiplicities of the infinite poles of `P(λ)`, returned in `ip`, are the positive differences between the multiplicities of
the infinite eigenvalues of `M-λN` and the degree of `P(λ)` [2].
The number of finite poles in `val` is equal to the number of finite eigenvalues of `M-λN` (returned in `KRInfo.nf`),
while the number of infinite poles in `val` is the sum of multiplicites in `ip`.
The multiplicities of the infinite eigenvalues of `Q(λ)` are returned in `id`.

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

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

[1] F. De Terán, F. M. Dopico, D. S. Mackey, Spectral equivalence of polynomial matrices and
the Index Sum Theorem, Linear Algebra and Its Applications, vol. 459, pp. 264-333, 2014.

[2] A. Varga, On computing the Kronecker structure of polynomial matrices using Julia, June 2020,
[arXiv:2006.06825](https://arxiv.org/pdf/2006.06825).
"""
function pmpoles(P::AbstractArray{T,3}; kwargs...) where T
   p, m, k1 = size(P)
   Ptilde = zeros(T,p+m,p+m,k1);
   Ptilde[1:p,1:m,:] = P;
   Ptilde[1:p,m+1:m+p,1] = Matrix{T}(I,p,p);
   Ptilde[p+1:p+m,1:m,1] = Matrix{T}(I,m,m);
   val, ip, kinfo = pmzeros(Ptilde; kwargs...)
   return val, ip, kinfo.id
end
pmpoles(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       pmpoles(poly2pm(P); kwargs...)
"""
    pmzeros1(P; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the finite and infinite zeros of the polynomial matrix `P(λ)` in `val`,
the multiplicities of infinite zeros in `iz` and information on the
Kronecker-structure of the underlying strongly irreducible pencil based linearization of `P(λ)`
in the `KRInfo` object.

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

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The computation of the zeros is performed by first building a strongly irreducible pencil based
linearization of `P(λ)` [1] of order `ν` as a structured system matrix pencil

              | A-λE | B-λF |
     M - λN = |------|------| ,
              | C-λG | D-λH |

such that `P(λ) = (C-λG)*inv(λE-A)*(B-λF)+D-λH`, and then reducing the pencil `M-λN` to an
appropriate Kronecker-like form (KLF) which exhibits information on its Kronecker structure.
The left and right Kronecker indices of `M-λN` and `P(λ)` are the same [2] and are
returned in `KRInfo.rki` and `KRInfo.rki`, respectively.
The multiplicities of the infinite zeros of `M-λN` and of `P(λ)` are the same [2] and are returned in `iz`.
The number of finite zeros in `val` is equal to the number of finite eigenvalues of `M-λN`
(returned in `KRInfo.nf`), while the number of infinite zeros in `val` is the sum of multiplicites in `iz`.
The normal rank of `P(λ)` can be evaluated as `KRInfo.nrank - ν`.

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

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

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

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.
"""
function pmzeros1(P::AbstractArray{T,3}; fast = false, atol::Real = zero(real(T)),
                  rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   if size(P,1) <= size(P,2)
      A, E, B, F, C, G, D, H, = lpsminreal(pm2lps(P; obs = true)...; fast = fast,
                                           atol1 = atol, atol2 = atol, rtol = rtol, obs = false)
   else
      A, E, B, F, C, G, D, H, = lpsminreal(pm2lps(P; contr = true)...; fast = fast,
                                           atol1 = atol, atol2 = atol, rtol = rtol, contr = false)
   end
   val, iz, info = pzeros([A B; C D], [E F; G H]; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   info.nrank -= size(A,1)
   return val, iz, info
end
pmzeros1(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) where {T} =
       pmzeros1(poly2pm(P); kwargs...)
"""
    pmpoles1(P; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, ip, id)

Return the infinite poles of the polynomial matrix `P(λ)` in `val`, the multiplicities of infinite poles in `ip` and
the infinite elementary divisors of the pole pencil `Sp(λ)` of the underlying strongly irreducible
pencil based linearization of  `Q(λ) := [P(λ) I; I 0]` in `id`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

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

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

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

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

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.
"""
function pmpoles1(P::AbstractArray{T,3}; fast = false, atol::Real = zero(real(T)),
                  rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   p = size(P,1)
   m = size(P,2)
   if p <= m
      A, E, B, F, C, G, D, H, = lpsminreal(pm2lps(P; obs = true)...; fast = fast,
                                           atol1 = atol, atol2 = atol, rtol = rtol, obs = false)
   else
      A, E, B, F, C, G, D, H, = lpsminreal(pm2lps(P; contr = true)...; fast = fast,
                                           atol1 = atol, atol2 = atol, rtol = rtol, contr = false)
   end
   n = size(A,1)
   M = [A zeros(T,n,p+m); zeros(T,p,n+m) I; zeros(T,m,n) I zeros(T,m,p)]
   N = [E F zeros(T,n,p); G H zeros(T,p,p); zeros(T,m,n+m+p)]
   val, iz, info = pzeros(M, N; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   return val, iz, info.id
end
pmpoles1(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       pmpoles1(poly2pm(P); kwargs...)
"""
      pmzeros2(P; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, iz, KRInfo)

Return the finite and infinite zeros of the polynomial matrix `P(λ)` in `val`,
the multiplicities of infinite zeros in `iz` and information on the
Kronecker-structure of the underlying irreducible descriptor system based linearization of `P(λ)`
in the `KRInfo` object.

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

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The computation of the zeros is performed by first building an irreducible descriptor system based
linearization of `P(λ)` [1] of order `ν` as a structured system matrix pencil

              | A-λE | B |
     M - λN = |------|---| ,
              |  C   | D |

such that `P(λ) = C*inv(λE-A)*B+D`, and then reducing the pencil `M-λN` to an
appropriate Kronecker-like form (KLF) which exhibits information on its Kronecker structure.
The left and right Kronecker indices of `M-λN` and `P(λ)` are the same [2] and are
returned in `KRInfo.rki` and `KRInfo.rki`, respectively.
The multiplicities of the infinite zeros of `M-λN` and of `P(λ)` are the same [2] and are returned in `iz`.
The number of finite zeros in `val` is equal to the number of finite eigenvalues of `M-λN`
(returned in `KRInfo.nf`), while the number of infinite zeros in `val` is the sum of multiplicites in `iz`.
The normal rank of `P(λ)` can be evaluated as `KRInfo.nrank - ν`.

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

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

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

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.
"""
function pmzeros2(P::AbstractArray{T,3}; fast = false, atol::Real = zero(real(T)),
                  rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   sys = pm2ls(P; contr = true, obs = true, fast = fast, atol = atol, rtol = rtol)
   val, iz, info = spzeros(sys...; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   info.nrank -= size(sys[1],1)
   return val, iz, info
end
pmzeros2(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       pmzeros2(poly2pm(P); kwargs...)
"""
    pmpoles2(P; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> (val, ip, id)

Return the infinite poles of the polynomial matrix `P(λ)` in `val`,
the multiplicities of infinite poles in `ip` and the infinite elementary divisors of the pole pencil
of the underlying irreducible descriptor system based linearization of `P(λ)` in `id`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

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

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

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

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

[2] G. Verghese, Comments on ‘Properties of the system matrix of a generalized state-space system’,
Int. J. Control, Vol.31(5) (1980) 1007–1009.
"""
function pmpoles2(P::AbstractArray{T,3}; fast = false, atol::Real = zero(real(T)),
                  rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   A, E, _ = pm2ls(P; contr = true, obs = true, fast = fast, atol = atol, rtol = rtol)
   val, ip, info = pzeros(A, E; fast = fast, atol1 = atol, atol2 = atol, rtol = rtol)
   return val, ip, info.id
end
pmpoles2(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       pmpoles2(poly2pm(P); kwargs...)
"""
    pmroots(P; fast = false, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> val

Compute the roots of the determinant of the regular polynomial matrix `P(λ)` in `val`.

`P(λ)` can be specified as a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`,
for which the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The roots of `det(P(λ))` are computed as the finite eigenvalues of a companion form linearization `M-λN`of `P(λ)`.
The finite eigenvalues are computed by reducing the pencil `M-λN` to an appropriate Kronecker-like form (KLF)
exhibiting the spliting of the infinite and finite eigenvalue structures of the pencil `M-λN`.
The reduction is performed using orthonal similarity transformations and involves rank decisions based on rank revealing QR-decompositions with column pivoting,
if `fast = true`, or, the more reliable, SVD-decompositions, if `fast = false`. For efficiency purposes, the reduction is only
partially performed, without accumulating the performed orthogonal transformations.

The keyword arguments `atol`  and `rtol` specify the absolute and relative tolerances for the nonzero
coefficients of `P(λ)`,  respectively.
The default relative tolerance is `n*ϵ`, where `n` is the size of the smallest dimension of `P(λ)`, and `ϵ` is the
machine epsilon of the element type of coefficients of `P(λ)`.
"""
function pmroots(P::AbstractArray{T,3}; fast = false, atol::Real = zero(real(T)),
                 rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   ispmregular(P,atol = atol, rtol = rtol) || error("The polynomial matrix is not regular")
   val, iz, _ = pmzeros(P, fast = fast, atol = atol, rtol = rtol)
   return val[1:end-sum(iz)]
end
pmroots(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
        pmroots(poly2pm(P); kwargs...)
"""
       pmrank(P; fastrank = true, atol = 0, rtol = atol > 0 ? 0 : n*ϵ)

Determine the normal rank of a polynomial matrix `P(λ)`.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`). The normal rank of `P(λ)`
is the number of linearly independent rows or columns.

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

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

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of `P(λ)`. The default relative tolerance is `n*ϵ`, where `n` is the size of the
smallest dimension of `P(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.
"""
function pmrank(P::AbstractArray{T,3}; fastrank::Bool = true, atol::Real = zero(real(T)),
                rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   if fastrank
       return rank(pmeval(P,exp(rand()*im)), atol = atol, rtol = rtol)
   else
       p, m, k = size(P)
       d = pmdeg(P)
       if p < m
          M, N = pm2lpCF2(P)
          r = prank(M, N, atol1 = atol, atol2 = atol, rtol = rtol)
          return d < 2 ? r : r-p*(d-1)
       else
          M, N = pm2lpCF1(P)
          r = prank(M, N, atol1 = atol, atol2 = atol, rtol = rtol)
          return d < 2 ? r : r-m*(d-1)
       end
   end
end
pmrank(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       pmrank(poly2pm(P); kwargs...)
"""
    ispmregular(P; fastrank = true, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> Bool

Test whether the polynomial matrix `P(λ)` is regular (i.e., `P(λ)` is square and ``{\\small\\det(P(λ)) \\not\\equiv 0}``).
The underlying computational procedure checks that the normal rank of the square `P(λ)` is equal to its order.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.
f
If `fastrank = true`, the rank is evaluated by counting how many singular values of `P(γ)` have magnitude greater
than `max(atol, rtol*σ₁)`, where `σ₁` is the largest singular value of `P(γ)` and `γ` is a randomly generated
complex value of magnitude equal to one.
If `fastrank = false`, first a linearization of `P(λ)` is built in a companion form `M-λN` of order `n`
and then its normal rank `nr` is determined. `P(λ)` is regular if `nr = n`.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of `P(λ)`. The default relative tolerance is `n*ϵ`, where `n` is the size of the
smallest dimension of `P(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.
"""
function ispmregular(P::AbstractArray{T,3}; fastrank = true, atol::Real = zero(real(T)),
                     rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   p, m, = size(P)
   m == p || (return false)
   if fastrank
      return rank(pmeval(P,exp(rand()*im)), atol = atol, rtol = rtol) == m
   else
      p <= m ? MN = pm2lpCF1(P) : MN = pm2lpCF2(P)
      return isregular(MN..., atol1 = atol, atol2 = atol, rtol = rtol)
   end
end
ispmregular(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
       ispmregular(poly2pm(P); kwargs...)
"""
    ispmunimodular(P; fastrank = true, atol::Real = 0, rtol::Real = atol>0 ? 0 : n*ϵ) -> Bool

Test whether the polynomial matrix `P(λ)` is unimodular (i.e., `P(λ)` is square, regular and `det(P(λ)) == constant`).

The underlying computational procedure checks that `P(λ)` is square and its first companion form linearization
is regular and has no finite eigenvalues.

`P(λ)` is a grade `k` polynomial matrix of the form `P(λ) = P_1 + λ P_2 + ... + λ**k P_(k+1)`, for which
the coefficient matrices `P_i`, `i = 1, ..., k+1`, are stored in the 3-dimensional matrix `P`,
where `P[:,:,i]` contains the `i`-th coefficient matrix `P_i` (multiplying `λ**(i-1)`).

`P(λ)` can also be specified as a matrix, vector or scalar of elements of the `Polynomial` type
provided by the [Polynomials](https://github.com/JuliaMath/Polynomials.jl) package.

The keyword arguments `atol` and `rtol`, specify, respectively, the absolute and relative tolerance for the
nonzero coefficients of `P(λ)`. The default relative tolerance is `n*ϵ`, where `n` is the size of the
smallest dimension of `P(λ)` and `ϵ` is the machine epsilon of the element type of its coefficients.
"""
function ispmunimodular(P::AbstractArray{T,3}; atol::Real = zero(real(T)),
                        rtol::Real = (min(size(P)...)*eps(real(float(one(T)))))*iszero(atol)) where T
   p, m, = size(P)
   m == p || (return false)
   return isunimodular(pm2lpCF1(P)..., atol1 = atol, atol2 = atol, rtol = rtol)
end
ispmunimodular(P::Union{AbstractVecOrMat{<:Polynomial},Polynomial,Number,AbstractVecOrMat{<:Number}}; kwargs...) =
               ispmunimodular(poly2pm(P); kwargs...)


text:
"A `Token` can only be held by one async process at one time. Use `Base.take!(token)` to claim the token, `Base.put!(token)` to give the token back."
struct Token
    c::Channel{Nothing}
    Token() = let
        c = Channel{Nothing}(1)
        push!(c, nothing)
        new(c)
    end
end

Base.take!(token::Token) = Base.take!(token.c)
Base.put!(token::Token) = Base.put!(token.c, nothing)
Base.isready(token::Token) = Base.isready(token.c)
Base.wait(token::Token) = Base.put!(token.c, Base.take!(token.c))

function withtoken(f::Function, token::Token)
    take!(token)
    result = try
        f()
    catch e
        rethrow(e)
    finally
        put!(token)
    end
    result
end

###

"Track whether some task needs to be done. `request!(requestqueue)` will make sure that the task is done at least once after calling it. Multiple calls might get bundled into one."
mutable struct RequestQueue
    is_processing::Bool
    c::Channel{Nothing}
    RequestQueue() = new(false, Channel{Nothing}(1))
end

"Give a function (with no arguments) that should be called after a request."
function process(f::Function, requestqueue::RequestQueue)
    @assert !requestqueue.is_processing
    requestqueue.is_processing = true
    while true
        take!(requestqueue.c)
        f()
    end
end

function request!(queue::RequestQueue)
    if isready(queue.c)
        push!(queue.c, nothing)
    end
end

###

mutable struct Promise{T}
    value::Union{Nothing,Some{T}}
    task::Union{Nothing,Task}
end

"
    Promise{T}(f::Function)

Run `f` asynchronously, and return a `Promise` to its result of type `T`. Call `fetch` on the returned `Promise` to await the result.

It's just like a `Task`, except the result is a type parameter.

# Example

```
julia> p = Promise() do
    sleep(5)
    1 + 2
end;

julia> fetch(p)
3
```

"
function Promise{T}(f::Function) where T
    p = Promise{T}(nothing, nothing)
    p.task = @async begin
        p.value = Some(f())
    end
    return p
end
Promise(f::Function) = Promise{Any}(f)

function Base.fetch(p::Promise{T})::T where T
	wait(p.task)
	something(p.value)
end

"Like @async except it prints errors to the terminal. 👶"
macro asynclog(expr)
	quote
		@async begin
			# because this is being run asynchronously, we need to catch exceptions manually
			try
				$(esc(expr))
			catch ex
				bt = stacktrace(catch_backtrace())
				showerror(stderr, ex, bt)
				rethrow(ex)
			end
		end
	end
end

ext: jl


text:
using DelimitedFiles
using PyPlot
using RAFF

"""

    draw_problem(M; raff_output=nothing, model_str="logistic",
                 datafile="/tmp/output.txt")

    draw_problem(datafile::String="/tmp/output.txt"; kwargs...)

Draw the problem data given by a (`m`x`3`) `M` matrix. By default it
is assumed that the model is given by the logistic model. If no
arguments are given it is assumed that the data file is given in file
`/tmp/output.txt`.

If a [RAFFOutput](@ref) object is provided, then it plots the model
found and also the outliers (true and false positives).

Optional arguments:

  - `raff_output`: [RAFFOutput](@ref) object with the solution
    obtained
  - `model_str`: a string with the name of the model to be used to
    plot the solution. See [model_list](@ref) for details.

"""
function draw_problem(M; raff_output=nothing, model_str="logistic",
                      θsol=nothing)

    x = M[:, 1]
    y = M[:, 2]
    co = M[:, 3]

    true_outliers = findall(co .!= 0.0)

    PyPlot.scatter(x[co .== 0.0], y[co .== 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="o", s=50.0, linewidths=0.2)

    PyPlot.scatter(x[co .!= 0.0], y[co .!= 0.0], color=PyPlot.cm."Pastel1"(2.0/9.0),
                   marker="^", s=50.0, linewidths=0.2, label="Outliers")

    if raff_output != nothing

        n, model, modelstr = RAFF.model_list[model_str]

        modl1 = (x) -> model(x, raff_output.solution)

        t = minimum(x):0.01:maximum(x)
        PyPlot.plot(t, modl1.(t), color=PyPlot.cm."Set1"(2.0/9.0))

        # Draw outliers found by RAFF

        true_positives = intersect(true_outliers, raff_output.outliers)
        false_positives = setdiff(raff_output.outliers, true_positives)

        PyPlot.scatter(x[false_positives], y[false_positives],
                       color=PyPlot.cm."Pastel1"(0.0/9.0), marker="o",
                       linewidths=0.2, edgecolors="k", s=50.0, label="False positives")

        PyPlot.scatter(x[true_positives], y[true_positives],
                       color=PyPlot.cm."Pastel1"(0.0/9.0), marker="^",
                       s=50.0, linewidths=0.2, edgecolors="k", label="Identified outliers")

    end

    # Plot the true solution, if available
    if θsol != nothing

        modl2 = (x) -> model(x, θsol)

        PyPlot.plot(t, modl2.(t), color=PyPlot.cm."Set1"(1.0/9.0),
                    linestyle="--")

    end

    PyPlot.legend(loc="best")

    PyPlot.show()

    PyPlot.savefig("/tmp/figure.png", DPI=150)

end

function draw_problem(datafile::String="/tmp/output.txt"; kwargs...)

    fp = open(datafile, "r")

    N = parse(Int, readline(fp))

    M = readdlm(fp)

    close(fp)

    draw_problem(M; kwargs...)

end


text:
function julia_url(version)
    os = @compat @static VERSION >= v"0.5.0" ? Base.Sys.KERNEL : Base.OS_NAME
    arch = @compat @static VERSION >= v"0.5.0" ? Base.Sys.WORD_SIZE : Base.WORD_SIZE
    julia_url(version, os, arch)
end

function julia_url(version::VersionNumber, os::Symbol, arch::Integer)
    # Cannibalized from https://github.com/travis-ci/travis-build/blob/master/lib/travis/build/script/julia.rb
    if os === :Linux && arch == 64
        os_arch = "linux/x64"
        ext = "linux-x86_64.tar.gz"
        nightly_ext = "linux64.tar.gz"
    elseif os === :Linux && arch == 32
        os_arch = "linux/x32"
        ext = "linux-i686.tar.gz"
        nightly_ext = "linux32.tar.gz"
    elseif os === :Darwin && arch == 64
        os_arch = "mac/x64"
        ext = "mac64.dmg"
        nightly_ext = ext
    elseif os === :Windows && arch == 64
        os_arch = "winnt/x64"
        ext = "win64.exe"
        nightly_ext = ext
    elseif os === :Windows && arch == 32
        os_arch = "winnt/x86"
        ext = "win32.exe"
        nightly_ext = ext
    else
        error(logger, "Julia does not support $arch-bit $os")
    end

    # Note: We could probably get specific revisions if we really wanted to.
    # eg. https://status.julialang.org/download/osx10.7+ redirects to
    # https://s3.amazonaws.com/julianightlies/bin/osx/x64/0.5/julia-0.5.0-2bb94d6f99-osx.dmg

    major_minor = "$(version.major).$(version.minor)"
    future_release = Base.nextpatch(NIGHTLY)  # Note: Nightly expected to be a prerelease
    pre = version.prerelease
    if version >= future_release
        throw(ArgumentError("The version $version exceeds the latest known nightly build $NIGHTLY"))
    elseif version >= NIGHTLY
        # https://julialangnightlies-s3.julialang.org/bin/linux/x64/julia-latest-linux64.tar.gz
        url = "julialangnightlies-s3.julialang.org/bin/$os_arch/julia-latest-$nightly_ext"
    elseif version.patch == 0 && (version == Base.upperbound(version) || (length(pre) > 0 && pre[1] == "latest"))
        # https://julialang-s3.julialang.org/bin/linux/x64/0.6/julia-0.6-latest-linux-x86_64.tar.gz
        url = "julialang-s3.julialang.org/bin/$os_arch/$major_minor/julia-$major_minor-latest-$ext"
    else
        # https://julialang-s3.julialang.org/bin/linux/x64/0.5/julia-0.5.2-linux-x86_64.tar.gz
        url = "julialang-s3.julialang.org/bin/$os_arch/$major_minor/julia-$version-$ext"
    end

    return "https://$url"
end

function julia_url(version::AbstractString, os::Symbol, arch::Integer)
    if version == "nightly"
        ver = NIGHTLY
    elseif version == "release"
        latest_release = VersionNumber(NIGHTLY.major, NIGHTLY.minor - 1, 0, (), ("",))
        ver = latest_release
    else
        ver = VersionNumber(version)
    end

    return julia_url(ver, os, arch)
end

function log_output(cmd::Cmd)
    try
        stdout = Pipe()
        stderr = Pipe()

        p = spawn(pipeline(cmd, stdin=DevNull, stdout=stdout, stderr=stderr))

        close(stdout.in)
        close(stderr.in)

        # Julia logs print to stderr
        for l in readlines(stderr)
            if startswith(l, "INFO: ")
                debug(logger, join(split(l)[2:end], ' '))
            elseif startswith(l, "WARNING: ")
                warn(logger, join(split(l)[2:end], ' '))
            else
                debug(logger, l)
            end
        end

        for l in readlines(stdout)
            debug(logger, l)
        end
    catch e
        error(logger, e)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 7cb1f900-375e-11eb-3cc9-cf8a8fb4fc04
file = readlines("day_3.txt")

# ╔═╡ 8f19f0e8-375e-11eb-38fa-b79ed15de216
md"Getting all tree positions on the original map"

# ╔═╡ a6db569a-375e-11eb-3943-051a1aefc042
world = Dict()

# ╔═╡ 7ed9ada8-375f-11eb-1194-93353bb23be0
for (index,line) in enumerate(file)
	trees = eachmatch(r"#",line)

	trees_pos = map(x-> x.offset, trees)

	for pos in trees_pos
		world[(index-1,pos-1)] = true
	end
end

# ╔═╡ 21106f36-375f-11eb-387e-7b11364cc055
world

# ╔═╡ 895f12cc-375f-11eb-3b70-a5c0accd57bc
md"extending the map to all positions"

# ╔═╡ 9818c068-375f-11eb-12a6-2b99d75fa222
begin
	tam_column = size(file)[1]
	tam_row = length(file[1])  #the \n ocupy one space

	has_tree(x::Int64,y::Int64) = haskey(world, ( x,y % tam_row))
end

# ╔═╡ d2c55762-3761-11eb-2e27-4b019a678050
tam_column

# ╔═╡ 5b953f4e-3762-11eb-1e0a-fbd3f82403d4
tam_row

# ╔═╡ ea80ee60-375d-11eb-2c3f-fd90c04e37c3
toboggan_pos = [(x,3*x) for x in 1:tam_column-1]

# ╔═╡ 6e1d3084-3769-11eb-1523-9d7fc6c7b1e3
[has_tree(x,y) for (x,y) in toboggan_pos]

# ╔═╡ 72a220ca-3763-11eb-3622-d33c420c6fb5
sum([has_tree(x,y) for (x,y) in toboggan_pos])

# ╔═╡ 24f92922-3768-11eb-018d-b582fb9346ad
md"# Question 2"

# ╔═╡ 3e9f32cc-3768-11eb-0561-7fac93e67081
slope(dx,dy) = sum([has_tree(x,Int64(dy/dx*x)) for x in 0:dx:tam_column-1])

# ╔═╡ 58b6ee84-3768-11eb-0cb0-5f322f2c4cc8
slope(1,3)

# ╔═╡ 7366bd04-3768-11eb-2366-a19e90bb92ca
slope(1,1)*slope(1,3)*slope(1,5)*slope(1,7)*slope(2,1)

# ╔═╡ 99bc5a36-3768-11eb-23e0-2394b9706da1
[slope(1,1),slope(1,3),slope(1,5),slope(1,7),slope(2,1)]

# ╔═╡ Cell order:
# ╠═7cb1f900-375e-11eb-3cc9-cf8a8fb4fc04
# ╟─8f19f0e8-375e-11eb-38fa-b79ed15de216
# ╠═a6db569a-375e-11eb-3943-051a1aefc042
# ╠═7ed9ada8-375f-11eb-1194-93353bb23be0
# ╠═21106f36-375f-11eb-387e-7b11364cc055
# ╠═895f12cc-375f-11eb-3b70-a5c0accd57bc
# ╠═9818c068-375f-11eb-12a6-2b99d75fa222
# ╠═d2c55762-3761-11eb-2e27-4b019a678050
# ╠═5b953f4e-3762-11eb-1e0a-fbd3f82403d4
# ╠═ea80ee60-375d-11eb-2c3f-fd90c04e37c3
# ╠═6e1d3084-3769-11eb-1523-9d7fc6c7b1e3
# ╠═72a220ca-3763-11eb-3622-d33c420c6fb5
# ╠═24f92922-3768-11eb-018d-b582fb9346ad
# ╠═3e9f32cc-3768-11eb-0561-7fac93e67081
# ╠═58b6ee84-3768-11eb-0cb0-5f322f2c4cc8
# ╠═7366bd04-3768-11eb-2366-a19e90bb92ca
# ╠═99bc5a36-3768-11eb-23e0-2394b9706da1


text:
using FunctionWrappersWrappers
using Test

@testset "FunctionWrappersWrappers.jl" begin

  fwplus = FunctionWrappersWrapper(+, (Tuple{Float64,Float64}, Tuple{Int,Int}), (Float64,Int));
  @test fwplus(4.0, 8.0) === 12.0
  @test fwplus(4, 8) === 12

  fwexp2 = FunctionWrappersWrapper(exp2, (Tuple{Float64}, Tuple{Float32}, Tuple{Int}), (Float64,Float32,Float64));
  @test fwexp2(4.0) === 16.0
  @test fwexp2(4f0) === 16f0
  @test fwexp2(4) === 16.0

end

max_stars_repo_name: chriselrod/FunctionWrappersWrappers

max_issues_repo_name: chriselrod/FunctionWrappersWrappers

max_forks_repo_name: chriselrod/FunctionWrappersWrappers


text:
module PkgUtils

using Pkg
using TOML
using Chain

export get_environments, get_projectfiles, get_pkg_version

# inspired from Pkg.jl/src/Pkg.jl
pkgdir(depot=Pkg.depots1()) = joinpath(depot, "packages")
const MANIFEST_USAGE = joinpath(Pkg.logdir(), "manifest_usage.toml")

function get_environments()::Array{Union{Nothing, String}, 1}
    try
        _mani = TOML.parsefile(MANIFEST_USAGE)
        _mani = collect(keys(_mani))
        mask = .~occursin.(pkgdir(), _mani) .& .~occursin.(tempdir(), _mani)
        return dirname.(_mani[mask])
    catch
        @warn "There is no environment now. Return nothing."
        return nothing
    end
end

# Rewrite from DrWatson.jl/src/project_setup.jl
function findprojectfile(dir::Union{Nothing, AbstractString}=pwd())
    # look for project file in current dir and parents
    home = homedir()
    while true
        # Base.project_names return ("JuliaProject.toml", "Project.toml")
        for proj in Base.project_names
            file = joinpath(dir, proj)
            Base.isfile_casesensitive(file) && return file  # return file instead of dir
        end
        # bail at home directory
        dir == home && break
        # bail at root directory
        old, dir = dir, dirname(dir)
        dir == old && break
    end
    @warn "Could not find find a project file by recursively checking " *
          "given `path` and its parents. Returning `nothing` instead."
    return nothing
end

function get_projectfiles()::Array{Union{Nothing, String}, 1}
    envpath = get_environments()
    return findprojectfile.(envpath)
end

# Rewrite from https://bkamins.github.io/julialang/2021/02/27/pkg_version.html
function get_pkg_version(name::AbstractString)
    try
        return @chain Pkg.dependencies() begin
            values
            [x for x in _ if x.name == name]
            only
            _.version
        end
    catch
        @warn "The standard library does not have a version number. Return nothing."
        return nothing
    end
end

end


text:
mutable struct AnalyticVI{T,O<:AVIOptimizer} <: VariationalInference{T}
    ϵ::T # Convergence criteria
    n_iter::Int # Number of steps performed
    stoch::Bool # Flag for stochastic optimization
    batchsize::Int # Size of mini-batches
    ρ::T # Scaling coeff. for stoch. opt.
    HyperParametersUpdated::Bool # Flag for updating kernel matrices
    vi_opt::O # Local optimizer for the variational parameters

    function AnalyticVI{T}(ϵ::T, optimiser, batchsize::Int, stoch::Bool) where {T}
        vi_opt = AVIOptimizer(optimiser)
        return new{T,typeof(vi_opt)}(ϵ, 0, stoch, batchsize, one(T), true, vi_opt)
    end
end

"""
    AnalyticVI(;ϵ::T=1e-5)

Variational Inference solver for conjugate or conditionally conjugate likelihoods (non-gaussian are made conjugate via augmentation)
All data is used at each iteration (use [`AnalyticSVI`](@ref) for updates using minibatches)

## Keywords arguments
- `ϵ::Real` : convergence criteria
"""
AnalyticVI

"""
    AnalyticSVI(nMinibatch::Int; ϵ::T=1e-5, optimiser=RobbinsMonro())

Stochastic Variational Inference solver for conjugate or conditionally conjugate likelihoods (non-gaussian are made conjugate via augmentation).
See [`AnalyticVI`](@ref) for reference

## Arguments

- `nMinibatch::Integer` : Number of samples per mini-batches

## Keywords arguments

- `ϵ::T` : convergence criteria
- `optimiser` : Optimiser used for the variational updates. Should be an Optimiser object from the [Flux.jl](https://github.com/FluxML/Flux.jl) library, see list here [Optimisers](https://fluxml.ai/Flux.jl/stable/training/optimisers/) and on [this list](https://github.com/theogf/AugmentedGaussianProcesses.jl/tree/master/src/inference/optimisers.jl). Default is `RobbinsMonro()` (ρ=(τ+iter)^-κ)
"""
AnalyticSVI

function AnalyticVI(; ϵ::T=1e-5) where {T<:Real}
    return AnalyticVI{T}(ϵ, Descent(1.0), 0, false)
end

function AnalyticSVI(
    nMinibatch::Integer; ϵ::T=1e-5, optimiser=RobbinsMonro()
) where {T<:Real}
    return AnalyticVI{T}(ϵ, optimiser, nMinibatch, true)
end

function Base.show(io::IO, inference::AnalyticVI)
    return print(
        io, "Analytic$(is_stochastic(inference) ? " Stochastic" : "") Variational Inference"
    )
end

### Generic method for variational updates using analytical formulas ###
# Single output version
@traitfn function variational_updates(
    m::TGP, state, y
) where {T,L,TGP<:AbstractGPModel{T,L,<:AnalyticVI};!IsMultiOutput{TGP}}
    local_vars = local_updates!(
        state.local_vars,
        likelihood(m),
        y,
        mean_f(m, state.kernel_matrices),
        var_f(m, state.kernel_matrices),
    )
    state = merge(state, (; local_vars))
    natural_gradient!.(
        m.f,
        ∇E_μ(m, y, local_vars),
        ∇E_Σ(m, y, local_vars),
        ρ(inference(m)),
        opt(inference(m)),
        Zviews(m),
        state.kernel_matrices,
        state.opt_state,
    )
    state = global_update!(m, state)
    return merge(state, (; local_vars))
end
# Multioutput version
@traitfn function variational_updates(
    m::TGP, state, y
) where {T,L,TGP<:AbstractGPModel{T,L,<:AnalyticVI};IsMultiOutput{TGP}}
    local_vars =
        local_updates!.(
            state.local_vars,
            likelihood(m),
            y,
            mean_f(m, state.kernel_matrices),
            var_f(m, state.kernel_matrices),
        ) # Compute the local updates given the expectations of f
    state = merge(state, (; local_vars))
    natural_gradient!.(
        m.f,
        ∇E_μ(m, y, state),
        ∇E_Σ(m, y, state),
        ρ(inference(m)),
        opt(inference(m)),
        Zviews(m),
        state.kernel_matrices,
        state.opt_state,
    ) # Compute the natural gradients of u given the weighted sum of the gradient of f
    state = global_update!(m, state) # Update η₁ and η₂
    return state
end

# Wrappers for tuple of 1 element,
# when multiple f are needed, these methods can be simply overloaded
function local_updates!(
    local_vars,
    l::AbstractLikelihood,
    y,
    μ::Tuple{<:AbstractVector{T}},
    diagΣ::Tuple{<:AbstractVector{T}},
) where {T}
    return local_updates!(local_vars, l, y, only(μ), only(diagΣ))
end

# Coordinate ascent updates on the natural parameters ##
function natural_gradient!(
    gp::Union{VarLatent{T},TVarLatent{T}},
    ∇E_μ::AbstractVector,
    ∇E_Σ::AbstractVector,
    ::Real,
    ::AVIOptimizer,
    X::AbstractVector,
    kernel_matrices,
    opt_state,
) where {T}
    K = kernel_matrices.K
    gp.post.η₁ .= ∇E_μ .+ K \ pr_mean(gp, X)
    gp.post.η₂ .= -Symmetric(Diagonal(∇E_Σ) + inv(K) / 2)
    return gp
end

# Computation of the natural gradient for the natural parameters
function natural_gradient!(
    gp::SparseVarLatent{T},
    ∇E_μ::AbstractVector,
    ∇E_Σ::AbstractVector,
    ρ::Real,
    ::AVIOptimizer,
    Z::AbstractVector,
    kernel_matrices,
    opt_state,
) where {T}
    K, κ = kernel_matrices.K, kernel_matrices.κ
    opt_state.∇η₁ .= ∇η₁(∇E_μ, ρ, κ, K, pr_mean(gp, Z), nat1(gp))
    opt_state.∇η₂ .= ∇η₂(∇E_Σ, ρ, κ, K, nat2(gp))
    return gp
end

# Gradient of on the first natural parameter η₁ = Σ⁻¹μ
function ∇η₁(
    ∇μ::AbstractVector{T},
    ρ::Real,
    κ::AbstractMatrix{T},
    K::Cholesky{T,Matrix{T}},
    μ₀::AbstractVector,
    η₁::AbstractVector{T},
) where {T<:Real}
    return transpose(κ) * (ρ * ∇μ) + (K \ μ₀) - η₁
end

# Gradient of on the second natural parameter η₂ = -1/2Σ⁻¹
function ∇η₂(
    θ::AbstractVector{T},
    ρ::Real,
    κ::AbstractMatrix{<:Real},
    K::Cholesky{T,Matrix{T}},
    η₂::Symmetric{T,Matrix{T}},
) where {T<:Real}
    return -(ρκdiagθκ(ρ, κ, θ) + inv(K) / 2) - η₂
end

# Natural gradient for the ONLINE model (OSVGP) #
function natural_gradient!(
    gp::OnlineVarLatent{T},
    ∇E_μ::AbstractVector{T},
    ∇E_Σ::AbstractVector{T},
    ::Real,
    ::AVIOptimizer,
    Z::AbstractVector,
    kernel_matrices,
    opt_state,
) where {T}
    K = kernel_matrices.K
    κ = kernel_matrices.κ
    κₐ = kernel_matrices.κₐ
    previous_gp = opt_state.previous_gp
    prevη₁ = previous_gp.prevη₁
    invDₐ = previous_gp.invDₐ
    gp.post.η₁ = K \ pr_mean(gp, Z) + transpose(κ) * ∇E_μ + transpose(κₐ) * prevη₁
    gp.post.η₂ =
        -Symmetric(ρκdiagθκ(1.0, κ, ∇E_Σ) + transpose(κₐ) * invDₐ * κₐ / 2 + inv(K) / 2)
    return gp
end

# Once the natural parameters have been updated we need to convert back to μ and Σ
@traitfn function global_update!(
    model::TGP, state
) where {T,L,TGP<:AbstractGPModel{T,L,<:AnalyticVI};IsFull{TGP}}
    global_update!.(model.f)
    return state
end

global_update!(gp::VarLatent, ::AVIOptimizer, ::AnalyticVI) = global_update!(gp)

function global_update!(model::OnlineSVGP, state)
    global_update!.(model.f)
    return state
end

#Update of the natural parameters and conversion from natural to standard distribution parameters
@traitfn function global_update!(
    model::TGP, state
) where {T,L,TGP<:AbstractGPModel{T,L,<:AnalyticVI};!IsFull{TGP}}
    opt_state =
        global_update!.(model.f, inference(model).vi_opt, inference(model), state.opt_state)
    return merge(state, (; opt_state))
end

function global_update!(gp::SparseVarLatent, opt::AVIOptimizer, i::AnalyticVI, opt_state)
    if is_stochastic(i)
        state_η₁, Δ₁ = Optimisers.apply(
            opt.optimiser, opt_state.state_η₁, nat1(gp), opt_state.∇η₁
        )
        state_η₂, Δ₂ = Optimisers.apply(
            opt.optimiser, opt_state.state_η₂, nat2(gp).data, opt_state.∇η₂
        )
        gp.post.η₁ .+= Δ₁
        gp.post.η₂ .= Symmetric(Δ₂) + nat2(gp)
        opt_state = merge(opt_state, (; state_η₁, state_η₂))
    else
        gp.post.η₁ .+= opt_state.∇η₁
        gp.post.η₂ .= Symmetric(opt_state.∇η₂ + nat2(gp))
    end
    global_update!(gp)
    return opt_state
end

function global_update!(gp::OnlineVarLatent, ::Any, ::Any)
    return global_update!(gp)
end

# Computation of the ELBO for all model
# There are 4 parts : the (augmented) log-likelihood, the Gaussian KL divergence
# the augmented variable KL divergence, some eventual additional part (like in the online case)
@traitfn function ELBO(
    model::TGP, state::NamedTuple, y
) where {T,L,TGP<:AbstractGPModel{T,L,<:AnalyticVI};!IsMultiOutput{TGP}}
    tot = zero(T)
    tot +=
        ρ(inference(model)) * expec_loglikelihood(
            likelihood(model),
            inference(model),
            y,
            mean_f(model, state.kernel_matrices),
            var_f(model, state.kernel_matrices),
            state.local_vars,
        )
    tot -= GaussianKL(model, state)
    tot -= Zygote.@ignore(
        ρ(inference(model)) * AugmentedKL(likelihood(model), state.local_vars, y)
    )
    tot -= extraKL(model, state)
    return tot
end

# Multi-output version
@traitfn function ELBO(
    model::TGP, state::NamedTuple, y
) where {T,L,TGP<:AbstractGPModel{T,L,<:AnalyticVI};IsMultiOutput{TGP}}
    tot = zero(T)
    tot +=
        ρ(inference(model)) * mapreduce(
            expec_loglikelihood,
            +,
            likelihood(model),
            inference(model),
            y,
            mean_f(model, state.kernel_matrices),
            var_f(model, state.kernel_matrices),
            state.local_vars,
        )
    tot -= GaussianKL(model, state)
    tot -= Zygote.@ignore(
        sum(ρ(inference(model)) * AugmentedKL.(likelihood(model), state.local_vars, y))
    )
    return tot
end

function expec_loglikelihood(
    l::AbstractLikelihood,
    i::AnalyticVI,
    y,
    μ::Tuple{<:AbstractVector{T}},
    diagΣ::Tuple{<:AbstractVector{T}},
    state,
) where {T}
    return expec_loglikelihood(l, i, y, only(μ), only(diagΣ), state)
end


text:
const USE_GPU = haskey(ENV, "USE_GPU") ? parse(Bool, ENV["USE_GPU"]) : false
const do_viz  = haskey(ENV, "DO_VIZ")  ? parse(Bool, ENV["DO_VIZ"])  : false
const do_save = haskey(ENV, "DO_SAVE") ? parse(Bool, ENV["DO_SAVE"]) : false
using ParallelRandomFields
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
    using ParallelRandomFields.grf2D_CUDA
else
    @init_parallel_stencil(Threads, Float64, 2)
    using ParallelRandomFields.grf2D_Threads
end

using ImplicitGlobalGrid, MAT, Plots
import MPI

@views inn(A) = A[2:end-1,2:end-1]

@views function generate_grf2D_multixpu()
    cov_typ = "expon"
    # Physics
    lx, ly  = 100.0, 100.0  # domain size
    sf      = 1.0           # standard deviation
    # -- exponential setup
    cl_e    = (10.0, 8.0)   # correlation lengths in [x,y]
    # -- gaussian setup
    cl_g    = 5.0           # correlation length isotropic
    k_m     = 100.0         # maximum value of the wave number
    # Numerics
    nh      = 10000         # inner parameter, number of harmonics
    nx, ny  = 32, 32        # numerical grid resolution
    me, dims, nprocs, coords, comm = init_global_grid(nx, ny, 1) # MPI initialisation
    @static if USE_GPU select_device() end    # select one GPU per MPI local rank (if >1 GPU per node)
    # Derived numerics
    dx, dy  = lx/nx_g(), ly/ny_g()    # cell sizes
    # Array allocation
    Yf      = @zeros(nx, ny)
    # Visu init
    if do_viz || do_save
        nx_v, ny_v = (nx-2)*dims[1], (ny-2)*dims[2]
        if (nx_v*ny_v*sizeof(Data.Number) > 0.8*Sys.free_memory()) error("Not enough memory for visualization.") end
        Yf_v     = zeros(nx_v, ny_v) # global array for visu
        Yf_inn   = zeros(nx-2, ny-2) # no halo local array for visu
        y_sl     = Int(ceil((ny_g()-2)/2))
        Xi_g, Yi_g = LinRange(dx+dx/2, lx-dx-dx/2, nx_v), LinRange(dy+dy/2, ly-dy-dy/2, ny_v) # inner points only

    end
    if cov_typ=="expon"
        # Generate the 2D exponential covariance function
        grf2D_expon!(Yf, sf, cl_e, nh, nx, ny, dx, dy; me=me, co1=coords[1], co2=coords[2], do_reset=true)
    elseif cov_typ=="gauss"
        # Generate the 2D Gaussian covariance function
        grf2D_gauss!(Yf, sf, cl_g, nh, k_m, nx, ny, dx, dy; me=me, co1=coords[1], co2=coords[2], do_reset=true)
    else
        error("Undefined covariance function")
    end
    # Visualisation
    if do_viz || do_save
        Yf_inn .= inn(Yf); gather!(Yf_inn, Yf_v)
        if do_viz && me==0
            heatmap(Xi_g, Yi_g, Yf_v', aspect_ratio=1, xlims=(Xi_g[1],Xi_g[end]), ylims=(Yi_g[1],Yi_g[end]), c=:hot, title="2D RandomField")
            savefig("grf2D_multixpu_$(nx_v)x$(ny_v).png")
        end
        if do_save && me==0
            file = matopen("grf2D_multixpu_$(cov_typ).mat", "w"); write(file, "grf2D", Yf_v); close(file)
        end
    end
    return
end

generate_grf2D_multixpu()


text:
#!/bin/sh
# Step 1: change/make directory
cd /mnt/d/GitHub/AEM7130_test
mkdir create_files
cd create_files

# Step 2: create sub-sample (first 5) of adult data
for i in $"head -n 5 ../adult.data";
do $i >> sub_adult_data.txt;
done

# Step 3: split sub_adult_data.txt by each row
split --numeric-suffixes=1 --additional-suffix=.txt -l 1 sub_adult_data.txt file

# Step 4: rename files
rename 's/file/file_/' *

# Step 5 append all files
for file in $(ls *.txt);
do cat $file >> new_data_set.csv;
done

# Step 6 count how many males, append the number to output.txt
grep Male new_data_set.csv | wc -l >> output.txt

# Step 7 count how many unique entries for column 7, append to output.txt
cut -f7 -d" " new_data_set.csv | sort -u | wc -l >> output.txt

# Step 7 remove all created files except for output.txt
find . ! -name output.txt -delete

# To execute bash file: chmod u+x hw1_q3.sh then bash hw1_q3.sh


text:
#=
include("move_suicidal.jl")     # difficulty level 1
include("move_random.jl")       # difficulty level 2
include("move_protracted.jl")   # difficulty level 3
include("move_normal.jl")       # difficulty level 4
include("move_hard.jl")         # difficulty level 5
=#
include("move_helpers.jl")

function updateBoard(oldBoard, movesTable)

    #duplicate the game board recursively using i-1 version and set the gametype depending on board size
    board = duplicateBoard(oldBoard)
    if length(board) == 5
        gameType = "mini"
    elseif length(board) == 9
        gameType = "standard"
    elseif length(board) == 12
        gameType = "chu"
    elseif length(board) == 16
        gameType = "ten"
    end

    # holds all the pieces that the even player (player 1) has captured
    capturedByEven = Array{piece,1}
    # holds all the pieces that the odd player (player 2) has captured
    capturedByOdd = Array{piece,1}

    #keep track of the current move, the ith move, and if the game is still being played
    currentMove = 1
    movesPlayed = length(movesTable[1])
    gameActive = true # flag for whether game has been resigned
    while gameActive && (currentMove <= movesPlayed)
        #displayBoard(board) # for debugging
        #=
        get(movesTable[col][row]) to retrieve move data
         - row number is identical to move number
         - col 1: move_number
         - col 2: move_type (move, resign or drop)
         - col 3: sourcex
         - col 4: sourcey
         - col 5: targetx
         - col 6: targety
         - col 7: option (! for move ending in promotion, piece name for drop move)
         - col 8: i_am_cheating (non-NULL if cheating)
         - col 9: targetx2
         - col 10: targety2
        =#

        side = currentMove % 2
        moveType = get(movesTable[2][currentMove])
        sourceX = movesTable[3][currentMove]
        sourceY = movesTable[4][currentMove]
        targetX = movesTable[5][currentMove]
        targetY = movesTable[6][currentMove]
        option = movesTable[7][currentMove]
        targetX2 = movesTable[9][currentMove]
        targetY2 = movesTable[10][currentMove]
        targetX3 = movesTable[11][currentMove]
        targetY3 = movesTable[12][currentMove]

        #conditionals for actions based on the variable moveType
        if moveType == "move"
            movePiece(gameType, board, currentMove, sourceX, sourceY, targetX, targetY, option, targetX2, targetY2, targetX3, targetY3)
        elseif moveType == "drop"
            dropPiece(board, currentMove, option, targetX, targetY)
        elseif moveType == "resign"
            gameActive = false
        end

        currentMove += 1
    end
    return board, capturedByEven, capturedByOdd, gameActive
end

function simulateMove(board1, move_number, sourceX, sourceY, targetX, targetY)
    board = duplicateBoard(board1)
    board[targetY][targetX] = board[sourceY][sourceX] # move piece at source to target, removing piece in process
    board[targetY][targetX].x = targetX # setting piece x
    board[targetY][targetX].y = targetY # setting piece y
    board[sourceY][sourceX] = piece("", "", 0, 0, -1, false) #set source space to empty

    deservesPromotion = promotion_check(board, sourceX, sourceY, targetX, targetY)
    if deservesPromotion
        board[targetY][targetX].name = uppercase(board[targetY][targetX].name)
    end
    return board, move_number
end

# player is either "sente" or "gote"; playerTime is the time to be enterd into the meta table
function updateMetaTime(gameFile, player, playerTime)
    SQLite.query(gameFile, "update meta set value='$playerTime' where key='$(player*"_time")';")
end

# inputs: SQLite database file and move object
# doesn't return anything
function makeMove(db, amove::move)
    move_number = amove.move_number
    move_type = amove.move_type
    sourceX = amove.sourcex
    sourceY = amove.sourcey
    targetX = amove.targetx
    targetY = amove.targety
    option = amove.option
    targetX2 = amove.targetx2
    targetY2 = amove.targety2
    targetX3 = amove.targetx3
    targetY3 = amove.targety3

    #db=SQLite.DB(gameFile)

    if move_type == "resign"
        SQLite.query(db, "insert into moves(move_number, move_type) values ($move_number, '$move_type')")
    elseif move_type == "drop"
        SQLite.query(db, "insert into moves(move_number, move_type, targetx, targety, option) values ($move_number, '$move_type',$targetX, $targetY, '$option')")
    elseif move_type == "move"
        # if single move
        if (targetX2 <= 0 || targetY2 <= 0) && (targetX3 <= 0 || targetY3 <= 0)
            SQLite.query(db, "insert into moves(move_number, move_type, sourcex, sourcey, targetx, targety, option) values ($move_number, '$move_type', $sourceX, $sourceY, $targetX, $targetY, '$option')")
        # if double move
        elseif targetX3 <= 0 || targetY3 <= 0
            SQLite.query(db, "insert into moves(move_number, move_type, sourcex, sourcey, targetx, targety, option, targetx2, targety2) values ($move_number, '$move_type', $sourceX, $sourceY, $targetX, $targetY, '$option', $targetX2, $targetY2)")
        # if triple move
        else
            SQLite.query(db, "insert into moves(move_number, move_type, sourcex, sourcey, targetx, targety, option, targetx2, targety2, targetx3, targety3) values ($move_number, '$move_type', $sourceX, $sourceY, $targetX, $targetY, '$option', $targetX2, $targetY2, $targetX3, $targetY3)")
        end
    end
end

function getUserMove(move_number)
    tic()
    # get players move
    moveMade = noMove
    moveMade.move_number = move_number
    while true
        println("\nPlease enter your move in the following format (omit any irrelevant elements)")
        println("\tmove_type sourceX sourceY targetX targetY targetX2 targetY2 targetX3 targetY3 option")
        print("Enter your move now: ")
        moveSelection = split(lowercase(chomp(readline(STDIN))))

        #println(moveSelection) # for debugging

        if length(moveSelection) == 1 && moveSelection[1] == "resign"
            moveMade.move_type = moveSelection[1]
            break
        elseif length(moveSelection) == 4 && moveSelection[1] == "drop"
            moveMade.move_type = moveSelection[1]
            moveMade.targetx = parse(Int8, moveSelection[2])
            moveMade.targety = parse(Int8, moveSelection[3])
            moveMade.option = moveSelection[4]
            break
        elseif length(moveSelection) > 4 && moveSelection[1] == "move"
            if length(moveSelection) == 5 # single move
                moveMade.move_type = moveSelection[1]
                moveMade.sourcex = parse(Int8, moveSelection[2])
                moveMade.sourcey = parse(Int8, moveSelection[3])
                moveMade.targetx = parse(Int8, moveSelection[4])
                moveMade.targety = parse(Int8, moveSelection[5])
                break
            elseif length(moveSelection) == 6 && moveSelection[6] == "!" # single move with promotion
                moveMade.move_type = moveSelection[1]
                moveMade.sourcex = parse(Int8, moveSelection[2])
                moveMade.sourcey = parse(Int8, moveSelection[3])
                moveMade.targetx = parse(Int8, moveSelection[4])
                moveMade.targety = parse(Int8, moveSelection[5])
                moveMade.option = moveSelection[6]
                break
            elseif length(moveSelection) == 7 # double move
                moveMade.move_type = moveSelection[1]
                moveMade.sourcex = parse(Int8, moveSelection[2])
                moveMade.sourcey = parse(Int8, moveSelection[3])
                moveMade.targetx = parse(Int8, moveSelection[4])
                moveMade.targety = parse(Int8, moveSelection[5])
                moveMade.targetx2 = parse(Int8, moveSelection[6])
                moveMade.targety2 = parse(Int8, moveSelection[7])
                break
            elseif length(moveSelection) == 8 && moveSelection[8] == "!"  # double move with promotion
                moveMade.move_type = moveSelection[1]
                moveMade.sourcex = parse(Int8, moveSelection[2])
                moveMade.sourcey = parse(Int8, moveSelection[3])
                moveMade.targetx = parse(Int8, moveSelection[4])
                moveMade.targety = parse(Int8, moveSelection[5])
                moveMade.targetx2 = parse(Int8, moveSelection[6])
                moveMade.targety2 = parse(Int8, moveSelection[7])
                moveMade.option = moveSelection[8]
                break
            elseif length(moveSelection) == 9 # triple move
                moveMade.move_type = moveSelection[1]
                moveMade.sourcex = parse(Int8, moveSelection[2])
                moveMade.sourcey = parse(Int8, moveSelection[3])
                moveMade.targetx = parse(Int8, moveSelection[4])
                moveMade.targety = parse(Int8, moveSelection[5])
                moveMade.targetx2 = parse(Int8, moveSelection[6])
                moveMade.targety2 = parse(Int8, moveSelection[7])
                moveMade.targetx3 = parse(Int8, moveSelection[8])
                moveMade.targety3 = parse(Int8, moveSelection[9])
                break
            elseif length(moveSelection) == 10 && moveSelection[10] == "!"  # triple move with promotion
                moveMade.move_type = moveSelection[1]
                moveMade.sourcex = parse(Int8, moveSelection[2])
                moveMade.sourcey = parse(Int8, moveSelection[3])
                moveMade.targetx = parse(Int8, moveSelection[4])
                moveMade.targety = parse(Int8, moveSelection[5])
                moveMade.targetx2 = parse(Int8, moveSelection[6])
                moveMade.targety2 = parse(Int8, moveSelection[7])
                moveMade.targetx3 = parse(Int8, moveSelection[8])
                moveMade.targety3 = parse(Int8, moveSelection[9])
                moveMade.option = moveSelection[10]
                break
            end
        end
        println("That is not a valid move selection. Please try again")
    end
    return moveMade, toq()
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Definition of constants.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorne, CA.
#
#   [2] http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/egm08_wgs84.html
#       Accessed on 2017-08-07.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# Julian Day of J2000.0 epoch.
const JD_J2000 = 2451545.0

# Earth Equatorial radius [m].
const R0 = 6378137.0

# Earth mean radius [m].
const Rm = 6371009.0

# Standard gravitational parameter for Earth [m^3/s^2].
const m0 = 3.986004418e14

# Perturbation terms based on EGM-08 standard gravitational model [1, pp. 1039].
const J2 = 1.08262617385222e-3
const J3 = 2.53241051856772e-6
const J4 = 1.61989759991697e-6

# Sun radius [m].
const Rs = 6.963e8

# Earth's orbit mean motion [rad/s]
const ne = (360.0/365.2421897)*pi/180/86400

# Earth's angular rotation [rad/s] without LOD correction.
const we =  7.292_115_146_706_979e-5

# Conversion factor from AU to m.
const au2m = 149597870700.0

# Sun radiation emitted [J/sec].
const sunRad = 3.826e26

# WGS-84 Data [2].
const a_wgs84  = 6378137.0
const f_wgs84  = 1/298.257223563
const b_wgs84  = a_wgs84*(1-f_wgs84)
const e_wgs84  = sqrt( (a_wgs84^2-b_wgs84^2)/a_wgs84^2 )
const el_wgs84 = sqrt( (a_wgs84^2-b_wgs84^2)/b_wgs84^2 )


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfCap::DataFrame, dfESRRev::DataFrame, dfResRevenue::DataFrame, dfChargingcost::DataFrame, dfPower::DataFrame, dfEnergyRevenue::DataFrame, dfSubRevenue::DataFrame, dfRegSubRevenue::DataFrame)

Function for writing net revenue of different generation technologies.
"""
function write_net_revenue(path::AbstractString, inputs::Dict, setup::Dict, EP::Model, dfCap::DataFrame, dfESRRev::DataFrame, dfResRevenue::DataFrame, dfChargingcost::DataFrame, dfPower::DataFrame, dfEnergyRevenue::DataFrame, dfSubRevenue::DataFrame, dfRegSubRevenue::DataFrame)
	dfGen = inputs["dfGen"]
	T = inputs["T"]     			# Number of time steps (hours)
	Z = inputs["Z"]     			# Number of zones
	G = inputs["G"]     			# Number of generators
	COMMIT = inputs["COMMIT"]		# Thermal units for unit commitment

	# Create a NetRevenue dataframe
 	dfNetRevenue = DataFrame(region = dfGen[!,:region], Resource = inputs["RESOURCES"], zone = dfGen[!,:Zone], Cluster = dfGen[!,:cluster], R_ID = dfGen[!,:R_ID])

	# Add investment cost to the dataframe
	dfNetRevenue.Inv_cost_MW = dfGen[!,:Inv_Cost_per_MWyr] .* dfCap[1:end-1,:NewCap]
	dfNetRevenue.Inv_cost_MWh = dfGen[!,:Inv_Cost_per_MWhyr] .* dfCap[1:end-1,:NewEnergyCap]
	if setup["ParameterScale"] == 1
		dfNetRevenue.Inv_cost_MWh = dfNetRevenue.Inv_cost_MWh * (ModelScalingFactor) # converting Million US$ to US$
		dfNetRevenue.Inv_cost_MW = dfNetRevenue.Inv_cost_MW * (ModelScalingFactor) # converting Million US$ to US$
	end

	# Add operations and maintenance cost to the dataframe
	dfNetRevenue.Fixed_OM_cost_MW = dfGen[!,:Fixed_OM_Cost_per_MWyr] .* dfCap[1:end-1,:EndCap]
 	dfNetRevenue.Fixed_OM_cost_MWh = dfGen[!,:Fixed_OM_Cost_per_MWhyr] .* dfCap[1:end-1,:EndEnergyCap]
 	dfNetRevenue.Var_OM_cost_out = (dfGen[!,:Var_OM_Cost_per_MWh]) .* dfPower[1:end-1,:AnnualSum]
	if setup["ParameterScale"] == 1
		dfNetRevenue.Fixed_OM_cost_MW = dfNetRevenue.Fixed_OM_cost_MW * (ModelScalingFactor) # converting Million US$ to US$
		dfNetRevenue.Fixed_OM_cost_MWh = dfNetRevenue.Fixed_OM_cost_MWh * (ModelScalingFactor) # converting Million US$ to US$
		dfNetRevenue.Var_OM_cost_out = dfNetRevenue.Var_OM_cost_out * (ModelScalingFactor) # converting Million US$ to US$
	end

	# Add fuel cost to the dataframe
	dfNetRevenue.Fuel_cost = zeros(size(dfNetRevenue, 1))
	for i in 1:G
		dfNetRevenue.Fuel_cost[i] = sum(inputs["C_Fuel_per_MWh"][i,:] .* inputs["omega"] .* value.(EP[:vP])[i,:])
	end
	if setup["ParameterScale"] == 1
		dfNetRevenue.Fuel_cost = dfNetRevenue.Fuel_cost * (ModelScalingFactor^2) # converting Million US$ to US$
	end

	# Add storage cost to the dataframe
	dfNetRevenue.Var_OM_cost_in = zeros(size(dfNetRevenue, 1))
 	for y in inputs["STOR_ALL"]
 		dfNetRevenue.Var_OM_cost_in[y] = dfGen[y,:Var_OM_Cost_per_MWh_In] * sum(inputs["omega"] .* value.(EP[:vCHARGE])[y,:])
 	end
	if setup["ParameterScale"] == 1
		dfNetRevenue.Var_OM_cost_in = dfNetRevenue.Var_OM_cost_in * (ModelScalingFactor^2) # converting Million US$ to US$
	end
	# Add start-up cost to the dataframe
	dfNetRevenue.StartCost = zeros(size(dfNetRevenue, 1))
 	if (setup["UCommit"]>=1)
 		for y in COMMIT #dfGen[!,:R_ID]
 			dfNetRevenue.StartCost[y] = sum(value.(EP[:eCStart])[y,:])
 		end
 	end
	if setup["ParameterScale"] == 1
		dfNetRevenue.StartCost = dfNetRevenue.StartCost * (ModelScalingFactor^2) # converting Million US$ to US$
	end
	# Add charge cost to the dataframe
	dfNetRevenue.Charge_cost = zeros(size(dfNetRevenue, 1))
	if has_duals(EP) == 1
		dfNetRevenue.Charge_cost = dfChargingcost[!,:AnnualSum] # Unit is confirmed to be US$
	end

	# Add energy and subsidy revenue to the dataframe
	dfNetRevenue.EnergyRevenue = zeros(size(dfNetRevenue, 1))
	dfNetRevenue.SubsidyRevenue = zeros(size(dfNetRevenue, 1))
	if has_duals(EP) == 1
		dfNetRevenue.EnergyRevenue = dfEnergyRevenue[!,:AnnualSum] # Unit is confirmed to be US$
	 	dfNetRevenue.SubsidyRevenue = dfSubRevenue[!,:SubsidyRevenue] # Unit is confirmed to be US$
	end

	# Add capacity revenue to the dataframe
	dfNetRevenue.ReserveMarginRevenue = zeros(size(dfNetRevenue, 1))
 	if setup["CapacityReserveMargin"] > 0 && has_duals(EP) == 1 # The unit is confirmed to be $
 		dfNetRevenue.ReserveMarginRevenue = dfResRevenue[!,:AnnualSum]
 	end

	# Add RPS/CES revenue to the dataframe
	dfNetRevenue.ESRRevenue = zeros(size(dfNetRevenue, 1))
 	if setup["EnergyShareRequirement"] > 0 && has_duals(EP) == 1 # The unit is confirmed to be $
 		dfNetRevenue.ESRRevenue = dfESRRev[!,:AnnualSum]
 	end

	# Calculate emissions cost
	dfNetRevenue.EmissionsCost = zeros(size(dfNetRevenue, 1))
 	if setup["CO2Cap"] >=1 && has_duals(EP) == 1
 		for y in 1:G
			dfNetRevenue.EmissionsCost[y] = 0.0
			for cap in 1:inputs["NCO2Cap"]
 				if dfGen[y,:Zone] in findall(x->x==1, inputs["dfCO2CapZones"][:,cap])
					if setup["CO2Cap"]==1 # Mass-based
						# Cost = sum(sum(emissions of gen y * dual(CO2 constraint[cap]) for z in Z) for cap in setup["NCO2"])
						dfNetRevenue.EmissionsCost[y] += sum(value(EP[:eEmissionsByPlant][y,t]) * inputs["omega"][t] * (-1) * dual(EP[:cCO2Emissions_systemwide][cap])
							for t in 1:T)

					elseif setup["CO2Cap"]==2 # Demand + Rate-based
						# Cost = sum(sum(emissions for zone z * dual(CO2 constraint[cap]) for z in Z) for cap in setup["NCO2"])
						dfNetRevenue.EmissionsCost[y] += sum(value(EP[:eEmissionsByPlant][y,t]) * inputs["omega"][t] * (-1) * dual(EP[:cCO2Emissions_systemwide][cap])
							for t in 1:T)
					elseif setup["CO2Cap"]==3 # Generation + Rate-based
						if y in union(inputs["THERM_ALL"],inputs["VRE"], inputs["VRE"],inputs["MUST_RUN"],inputs["HYDRO_RES"])
							# Cost = sum( sum(emissions - generatio for zone z * MaxCO2Rate for zone z for z in Z) * dual(CO2 constraint[cap] for cap in setup["NCO2"])
							dfNetRevenue.EmissionsCost[y] += sum( (value(EP[:eEmissionsByPlant][y,t]) - (value(EP[:vP][y,t]) * inputs["dfMaxCO2Rate"][z,cap])) * inputs["omega"][t] * (-1) * dual(EP[:cCO2Emissions_systemwide][cap])
								for z=dfGen[y,:Zone], t in 1:T)
						end
					end
				end
			end
 		end
		if setup["ParameterScale"] == 1
			dfNetRevenue[!,:EmissionsCost] = dfNetRevenue[!,:EmissionsCost] * (ModelScalingFactor^2) # converting Million US$ to US$
		end
 	end

	# Add regional technology subsidy revenue to the dataframe
	dfNetRevenue.RegSubsidyRevenue = zeros(size(dfNetRevenue, 1))
	if setup["MinCapReq"] >= 1 && has_duals(EP) == 1 # The unit is confirmed to be US$
		dfNetRevenue.RegSubsidyRevenue = dfRegSubRevenue[!,:SubsidyRevenue]
	end

	dfNetRevenue.Revenue =	dfNetRevenue.EnergyRevenue + dfNetRevenue.SubsidyRevenue + dfNetRevenue.ReserveMarginRevenue + dfNetRevenue.ESRRevenue + dfNetRevenue.RegSubsidyRevenue
	dfNetRevenue.Cost = dfNetRevenue.Inv_cost_MW + dfNetRevenue.Inv_cost_MWh + dfNetRevenue.Fixed_OM_cost_MW + dfNetRevenue.Fixed_OM_cost_MWh + dfNetRevenue.Var_OM_cost_out + dfNetRevenue.Var_OM_cost_in + dfNetRevenue.Fuel_cost + dfNetRevenue.Charge_cost + dfNetRevenue.EmissionsCost + dfNetRevenue.StartCost
	#dfNetRevenue.Cost = dfNetRevenue.Inv_cost_MW + dfNetRevenue.Inv_cost_MWh + dfNetRevenue.Fixed_OM_cost_MW + dfNetRevenue.Fixed_OM_cost_MWh + dfNetRevenue.Var_OM_cost_out + dfNetRevenue.Var_OM_cost_in + dfNetRevenue.Fuel_cost + dfNetRevenue.Charge_cost + dfNetRevenue.EmissionsCost + dfNetRevenue.StartCost
	dfNetRevenue.Profit = 	dfNetRevenue.Revenue - dfNetRevenue.Cost

	CSV.write(joinpath(path, "NetRevenue.csv"), dfNetRevenue)
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
module testUniverse

using TimeData
using DataFrames
using Distributions
using NumericExtensions
using Winston

include(string(Pkg.dir("AssetMgmt"), "/src/AssetMgmt.jl"))

println("\n Running bootstrap analysis\n")

###############
## real data ##
###############

fullData = readTimedata("/home/chris/research/asset_mgmt/data/datesLogRet.csv")
data = fullData[:, 1:end]

## dataFile = string(Pkg.dir("AssetMgmt"), "/data/discr_ret.csv")
## data = readTimedata(dataFile)

univ = AssetMgmt.Universe(data)

## get "true" underlying distribution
pdMat = PDMat(array(univ.covMatr))
d = MvNormal(array(univ.mus)[:], pdMat)

muAddon = 0.5
nMuGrid = 40
muGrid = AssetMgmt.getMuGrid(univ, muAddon, nMuGrid)
effPfs = AssetMgmt.effPf(univ, muGrid)
effMus = AssetMgmt.getPMean(effPfs, univ)
effVars = AssetMgmt.getPVar(effPfs, univ)

## plot true efficient frontier
p = FramedPlot(title="Bootstrapped efficient frontiers",
         xlabel="\\sigma",
         ylabel="\\mu")

add(p, Curve(sqrt(effVars), effMus[:], color="red"))
add(p, Points(sqrt(diag(array(univ.covMatr))), array(univ.mus)))

nBootstr = 250
nObs = 4600

realizedMus = Array(Float64, length(muGrid), nBootstr)
realizedVars = Array(Float64, length(muGrid), nBootstr)

tic();

@profile begin
for ii=1:nBootstr
## for ii=1:5
    simVals = rand(d, nObs)'

    mus = mean(simVals, 1)
    covMatr = cov(simVals)

    fakeUniv = AssetMgmt.Universe(mus, covMatr, names(univ))
    fakeMuGrid = AssetMgmt.getMuGrid(fakeUniv, muAddon, nMuGrid)
    fakeEffPfs = AssetMgmt.effPf(fakeUniv, fakeMuGrid)

    realizedMus[:, ii] = AssetMgmt.getPMean(fakeEffPfs, univ)
    realizedVars[:, ii] = AssetMgmt.getPVar(fakeEffPfs, univ)

    ## realizedSigmas = sqrt(realizedVars)
    add(p, Curve(sqrt(realizedVars[:, ii]), realizedMus[:, ii], color="blue"))
end
end

toc()

p

#######################################
## visualized bootstrapped frontiers ##
#######################################


text:
# LSTM model - input * weight, concatenated weights
function lstm(weight, bias, hidden, cell, input)
    gates   = hcat(input,hidden) * weight .+ bias
    hsize   = size(hidden,2)
    forget  = sigm(gates[:,1:hsize])
    ingate  = sigm(gates[:,1+hsize:2hsize])
    outgate = sigm(gates[:,1+2hsize:3hsize])
    change  = tanh(gates[:,1+3hsize:end])
    cell    = cell .* forget + ingate .* change
    hidden  = outgate .* tanh(cell)
    return (hidden,cell)
end

function logprob(output, ypred, mask=nothing)
    nrows,ncols = size(ypred)
    index = similar(output)
    @inbounds for i=1:length(output)
        index[i] = i + (output[i]-1)*nrows
    end
    o1 = logp(ypred,2)
    index = mask == nothing ? index : index[mask]
    o2 = o1[index]
    o3 = sum(o2)
    return o3
end

# attention mechanism
function att(w,a,h, o=Dict())
    # get sizes
    B,L,D  = size(a); H = size(h,2)
    attdrop = get(o, :attdrop, 0.0)

    # flatten region features
    aflat = reshape(a, B*L, D)

    # visual features projection
    aatt = aflat*w["wce"] # size(wce): (D,D), projection
    aatt = reshape(aatt, B, L, D)

    # hidden states projection
    hatt = h * w["whe"]
    hatt = reshape(hatt, B, 1, D) # size(whe): (H,D)

    # et calculation
    et = tanh(aatt .+ hatt)
    et = reshape(et, B*L, D)
    et = dropout(et, attdrop) * w["watt"] # watt: D,1
    et = reshape(et, B, L)

    # alpha = logp(et, 2)
    et1 = maximum(et, 2)
    et2 = et .- et1
    et3 = exp(et2)
    et4 = sum(et3, 2)
    et5 = et3 ./ et4

    alpha = exp(et5)
    alpha = alpha ./ sum(alpha,2)
    alpha = reshape(alpha, size(alpha)..., 1)

    context = alpha .* a # B,L,1 * B,L,D
    context = sum(context,2)
    context = reshape(context, size(context,1), size(context,3))
end

max_stars_repo_name: ilkerkesen/show-attend-and-tell

max_issues_repo_name: ilkerkesen/show-attend-and-tell

max_forks_repo_name: ilkerkesen/show-attend-and-tell


text:
push!(LOAD_PATH, homedir()*"/Projects/QuIPS/src")

using QuIPS

# ψ₀ = Ĥ₁ * [0,1] ⊗ [1,0] ⊗ [1,0]

quip = [
    # entangeling alice and bob's qubits
    (:H, 2),
    (:CX, (2, 3)),

    # preparing qubit to be teleported by alice to bob
    (:X, 1),
    (:H, 1),

    # teleportation algorithm
    (:CX, (1, 2)),
    (:H, 1),
    (:M, 1),
    (:M, 2)
]

QC = QCircuit(quip, 3)

QC |> run!

m₁ = Int(QC.out[1])
m₂ = Int(QC.out[2])

X̂(i) = Gate((:X, i))
Ẑ(i) = Gate((:Z, i))

operate!(QC, X̂(3)^m₂)
X̂′ = QC.ops[end]
print("  "); show(X̂′)

operate!(QC, Ẑ(3)^m₁)
Ẑ′ = QC.ops[end]
print("  "); show(Ẑ′)

println("\nQC.wfn:\n")
for i in QC.wfn println("  ", i) end

X = Gates().Us[:X]
H = Gates().Us[:H]

ψ′ = (X^m₁ * [1,0]) ⊗ (X^m₂ * [1,0]) ⊗ (H * X * [1,0])

const C = Complex{Float32}

println("\nψ′:\n")
for i in C.(ψ′) println("  ", i) end

M̂₃ = Measurement(3)

operate!(QC, M̂₃)

println("\nψ′′:\n")
for i in QC.wfn println("  ", i) end

println("\nm₃ = ", QC.out[3], "\n")

max_stars_repo_name: aarontrowbridge/QuIPS

max_issues_repo_name: aarontrowbridge/QuIPS

max_forks_repo_name: aarontrowbridge/QuIPS


text:
import Pkg

if Base.VERSION >= VersionNumber("1.1")
    Pkg.Registry.add(
        Pkg.RegistrySpec(
            name="BCBIRegistry",
            url="https://github.com/bcbi/BCBIRegistry.git",
            uuid="26a550a3-39fe-4af4-af6d-e8814c2b6dd9",)
        )
    Pkg.Registry.add("General")
    Pkg.Registry.update()
end

root = joinpath(dirname(@__FILE__))
@info("Updating required environment at $(root)")
Pkg.activate(root)
Pkg.instantiate()
Pkg.update()
if Base.VERSION >= VersionNumber("1.3")
    Pkg.precompile()
else
    Pkg.API.precompile()
end

required_environment_list = [
    ["environments", "datasets",],
    ["environments", "interop",],
    ["environments", "interop-seaborn",],
    ["environments", "plotting",],
    ["environments", "predictmd",],
    ]

optional_environment_list = [
    ["environments", "vegalite",],
    ]

for environment in required_environment_list
    location = joinpath(dirname(@__FILE__), environment...)
    @info("Updating required environment at $(location)")
    Pkg.activate(location)
    Pkg.instantiate()
    Pkg.update()
    if Base.VERSION >= VersionNumber("1.3")
        Pkg.precompile()
    else
        Pkg.API.precompile()
    end
end

for environment in optional_environment_list
    location = joinpath(dirname(@__FILE__), environment...)
    @info("Updating optional environment at $(location)")
    Pkg.activate(location)
    try
        Pkg.instantiate()
    catch ex
        @error("Encountered a non-fatal error", exception=ex, environment)
    end
    try
        Pkg.update()
    catch ex
        @error("Encountered a non-fatal error", exception=ex, environment)
    end
    if Base.VERSION >= VersionNumber("1.3")
        try
            Pkg.precompile()
        catch ex
            @error("Encountered a non-fatal error", exception=ex, environment)
        end
    else
        try
            Pkg.API.precompile()
        catch ex
            @error("Encountered a non-fatal error", exception=ex, environment)
        end
    end
end

max_stars_repo_name: bcbi/BCBI-base

max_issues_repo_name: bcbi/BCBI-base

max_forks_repo_name: bcbi/BCBI-base