text:
# Copyright (c) 2021, Eric Sabo
# All rights reserved.
#
# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

ENV["NEMO_PRINT_BANNER"] = "false"

using AbstractAlgebra
using Nemo

import Base: show, length, in, ⊆, /, *
import AbstractAlgebra: quo, VectorSpace

include("utils.jl")

abstract type AbstractCode end
abstract type AbstractLinearCode <: AbstractCode end

# can't get around including this here and not in weight_dist without modules
struct WeightEnumerator
    polynomial::Vector{Vector{Int64}}
    Type::String
end

mutable struct LinearCode <: AbstractLinearCode
    F::Union{FqNmodFiniteField, Nemo.GaloisField, AbstractAlgebra.GFField{Int64}}
    n::Integer
    k::Integer
    d::Union{Integer, Missing}
    G::Union{gfp_mat, fq_nmod_mat}
    Gorig::Union{gfp_mat, fq_nmod_mat, Missing}
    H::Union{gfp_mat, fq_nmod_mat}
    Horig::Union{gfp_mat, fq_nmod_mat, Missing}
    Gstand::Union{gfp_mat, fq_nmod_mat}
    Hstand::Union{gfp_mat, fq_nmod_mat}
    weightenum::Union{WeightEnumerator, Missing}
end

# need to make sure this isn't doing column swaps
function _standardform(G::Union{gfp_mat, fq_nmod_mat})
    # hnf(C.G) shoud also work in most situations
    _, Gstand = rref(G) # breaks on this line
    A = Gstand[:, (size(Gstand, 1) + 1):size(Gstand, 2)]
    T = MatrixSpace(base_ring(Gstand), size(Gstand, 2) - size(Gstand, 1), size(Gstand, 2) - size(Gstand, 1))
    Hstand = hcat(-transpose(A), T(1))
    return Gstand, Hstand
end

"""
    LinearCode(G::fq_nmod_mat, parity::Bool=false, verify::Bool=true)

Return a LinearCode object using the matrix `G` as the generator matrix.

If `G` is not full rank, a row-reduced form is computed for the generator matrix.
The dimension of the code is the number of rows of the full-rank matrix, and the
length the number of columns. If the optional paramater `parity` is set to `true`,
a LinearCode is built with `G` as the parity check matrix. If the optional parameter
`verify` is set to `true`, basic checks, such as checking `GH == 0` and checking
for column swaps in the standard form, will be done to ensure correctness.

# Arguments
* `G`: A matrix over a finite field of type `FqNmodFiniteField`.
* `parity`: A boolean
* `verify`: A boolean

# Notes
* At the moment, no convention is used for G = 0 and an error is thrown.
* Zero columns are not removed.
* Row reduction is based on Nemo's rref function. It should be monitored to make
  sure that this does not introduce column swapping in a future version.

# Examples
```julia
julia> F, _ = FiniteField(2, 1, "α");
julia> G = matrix(F, 4, 7, [1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1]);
julia> C = LinearCode(G)
[7, 4]_2 linear code.
Generator matrix: 4 × 7
        1 0 0 0 0 1 1
        0 1 0 0 1 0 1
        0 0 1 0 1 1 0
        0 0 0 1 1 1 1
```
"""
function LinearCode(G::Union{gfp_mat, fq_nmod_mat}, parity::Bool=false, verify::Bool=true)
    # if passed a VS, just ask for the basis
    !iszero(G) || error("Zero matrix passed into LinearCode constructor.")
    if typeof(base_ring(G)) != FqNmodFiniteField && typeof(base_ring(G)) != Nemo.GaloisField && typeof(base_ring(G)) != AbstractAlgebra.GFField{Int64}
        error("Codes over general rings are currently not supported, only fields.")
    end

    Gorig = G
    rk = rank(G)
    !iszero(rk) || error("Rank zero matrix passed into LinearCode constructor.")
    if rk < size(G, 1)
        _, G = rref(G)
    end

    del = Vector{Int64}()
    for r in 1:size(G, 1)
        if iszero(G[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        G = G[setdiff(1:size(G, 1), del), :]
    end

    # note the H here is transpose of the standard definition
    _, H = right_kernel(G)

    Gstand, Hstand = _standardform(G)
    if verify
        size(H', 1) == size(G, 2) - size(G, 1) || error("Parity check matrix is not size n - k.")
        !(!iszero(G * H) || !iszero(H' * G')) || error("Generator and parity check matrices are not transpose orthogonal.")
        for r in 1:size(Gstand, 1)
            iszero(Gstand[r, :] * H) || error("Column swap appeared in _standardform.")
        end
    end

    if !parity
        return LinearCode(base_ring(G), size(G, 2), size(G, 1), missing, G, Gorig,
            deepcopy(H'), missing, Gstand, Hstand, missing)
    else
        return LinearCode(base_ring(G), size(G, 2), size(G, 1), missing, deepcopy(H'), missing,
            G, Gorig, Hstand, Gstand, missing)
    end
end

# need to unify notation here to call minimumdistance but since this was removed
# from this file and put into weight_dist, keep this until exporting
function show(io::IO, C::AbstractLinearCode)
    if get(io, :compact, false)
        if ismissing(C.d)
            println(io, "[$(length(C)), $(dimension(C))]_$(order(field(C))) linear code.")
        else
            println(io, "[$(length(C)), $(dimension(C)), $(C.d)]_$(order(field(C))) linear code.")
        end
    else
        if ismissing(C.d)
            println(io, "[$(length(C)), $(dimension(C))]_$(order(field(C))) linear code.")
        else
            println(io, "[$(length(C)), $(dimension(C)), $(C.d)]_$(order(field(C))) linear code.")
        end
        println(io, "Generator matrix: $(dimension(C)) × $(length(C))")
        for i in 1:dimension(C)
            print(io, "\t")
            for j in 1:length(C)
                if j != length(C)
                    print(io, "$(C.G[i, j]) ")
                elseif j == length(C) && i != dimension(C)
                    println(io, "$(C.G[i, j])")
                else
                    print(io, "$(C.G[i, j])")
                end
            end
        end
    end
end

"""
    field(C::AbstractLinearCode)

Return the base ring of the generator matrix as a Nemo object.
"""
field(C::AbstractCode) = C.F

"""
    length(C::AbstractLinearCode)

Return the length of the code.
"""
length(C::AbstractCode) = C.n
# n(C::AbstracCode) = C.n

"""
    dimension(C::AbstractLinearCode)

Return the dimension of the linear code.
"""
dimension(C::AbstractLinearCode) = C.k
# k(C::AbstractLinearCode) = C.k

"""
    cardinality(C::AbstractLinearCode)

Return the cardinality of the linear code.

No size checking is done on the parameters of the code, returns a BitInt by default.
"""
cardinality(C::AbstractLinearCode) = BigInt(order(field(C)))^dimension(C)

"""
    rate(C::AbstractLinearCode)

Return the rate, `R = k/n', of the linear code.
"""
rate(C::AbstractLinearCode) = dimension(C) / length(C)

# """
#     minimumdistance(C::AbstractLinearCode)
#
# Return the minimum distance of the linear code if known, otherwise returns missing.
# """
# minimumdistance(C::AbstractLinearCode) = C.d

"""
    setminimumdistance(C::AbstractLinearCode, d::Integer)

Set the minimum distance of the linear code to `d`.

The only check done on the value of `d` is that `1 ≤ d ≤ n`.
"""
function setminimumdistance!(C::AbstractLinearCode, d::Integer)
    d > 0 && d <= length(C) || error("The minimum distance of a code must be ≥ 1; received: d = $d.")
    C.d = d
end

"""
    relativedistance(C::AbstractLinearCode)

Return the relative minimum distance, `δ = d / n` of the linear code if `d` is known,
otherwise errors.
"""
function relativedistance(C::AbstractLinearCode)
    !ismissing(C.d) || error("Unknown minimum distance for this code.")
    return C.d / C.n
end

"""
    generatormatrix(C::AbstractLinearCode, standform::Bool=false)

Return the generator matrix of the linear code.

If the optional parameter `standform` is set to `true`, the standard form of the
generator matrix is returned instead.
"""
function generatormatrix(C::AbstractLinearCode, standform::Bool=false)
    !standform || return C.Gstand
    return C.G
end

"""
    originalgeneratormatrix(C::AbstractLinearCode)

Return the original generator matrix used to create the linear code.

This should not be interpreted as any kind of generator to the code. For example,
if `G` is not a full-rank matrix, then this function acting on `C = LinearCode(G)`
will return `G` instead of the generator matrix for `C`. This may be used to diagnose
or understand the resulting code after its creation. An important use of this
function is in connection with standard code modification procedures such as
extending, puncturing, shortening, etc where this returns the matrix on which these
procedures have been applied.
"""
function originalgeneratormatrix(C::AbstractLinearCode)
    !ismissing(C.Gorig) || return C.G # is ths the behavior I want?
    return C.Gorig
end

"""
    paritycheckmatrix(C::AbstractLinearCode, standform::Bool=false)

Return the parity check matrix of the linear code.

If the optional parameter `standform` is set to `true`, the standard form of the
parity check matrix is returned instead.
"""
function paritycheckmatrix(C::AbstractLinearCode, standform::Bool=false)
    !standform || return C.Hstand
    return C.H
end

"""
    originalparitycheckmatrix(C::AbstractLinearCode)

Return the original parity check matrix used to create the linear code.

See `originalgeneratormatrix`.
"""
function originalparitycheckmatrix(C::AbstractLinearCode)
    !ismissing(C.Horig) || return C.H # is ths the behavior I want?
    return C.Horig
end

"""
    genus(C::AbstractLinearCode)

Return the genus, `n + 1 - k - d` the linear code.
"""
function genus(C::AbstractLinearCode)
    !ismissing(C.d) || error("Unknown minimum distance for this code.")
    return length(C) + 1 - dimension(C) - minimumdistance(C)
end

function Singletonbound(n::Integer, a::Integer)
    # d ≤ n - k + 1 or k ≤ n - d + 1
    if n >= 0 && a >= 0 && n >= a
        return n - a + 1
    else
        error("Invalid parameters for the Singleton bound. Received n = $n, k/d = $a")
    end
end

function Singletonbound(C::AbstractLinearCode)
    return Singletonbound(length(C), dimension(C))
end

"""
    genus(C::AbstractLinearCode)

Return whether or not the linear code is maximum distance separable (MDS).

A linear code is MDS if it saturates the Singleton bound, `d = n - k + 1`.
"""
function isMDS(C::AbstractLinearCode)
    minimumdistance(C) != Singletonbound(C) || return true
    return false
end

"""
    numbercorrectableerrors(C::AbstractLinearCode)

Return the number of correctable errors for the linear code.

The number of correctable errors is `t = floor((d - 1) / 2)`.
"""
function numbercorrectableerrors(C::AbstractLinearCode)
    !ismissing(minimumdistance(C)) || return -1 # what do I want to return here?
    return Int64(floor((minimumdistance(C) - 1) / 2))
end

"""
    encode(v::Union{fq_nmod_mat, Vector{Int64}}, C::AbstractLinearCode)

Return `v * G`, where `G` is the generator matrix of `C`.

# Arguments
* `v`: Either a `1 × k` or a `k × 1` vector.
"""
function encode(v::Union{gfp_mat, fq_nmod_mat}, C::AbstractLinearCode)
    !(size(v) != (1, dimension(C)) && size(v) != (dimension(C), 1)) ||
        error("Vector to be encoded is of incorrect dimension; expected length $(C.k), received: $(size(v)).")
    base_ring(v) == field(C) || error("Vector must have the same base ring as the generator matrix.")

    size(v, 1) != 1 || return v * generatormatrix(C)
    return v' * generatormatrix(C)
end

function encode(v::Vector{Int64}, C::AbstractLinearCode)
    length(v) == dimension(C) ||
        error("Vector to be encoded is of incorrect length; expected length $(C.k), received: $(size(v)).")
    return encode(matrix(field(C), v'), C)
end

"""
    syndrome(v::Union{fq_nmod_mat, Vector{Int64}}, C::AbstractLinearCode)

Return `H * v`, where `H` is the parity check matrix of `C`.

# Arguments
* `v`: Either a `1 × k` or a `k × 1` vector.
"""
function syndrome(v::Union{gfp_mat, fq_nmod_mat}, C::AbstractLinearCode)
    !(size(v) != (length(C), 1) && size(v) != (1, length(C))) ||
        error("Vector to be tested is of incorrect dimension; expected length $(C.n), received: $(size(v)).")
    base_ring(v) == field(C) || error("Vector must have the same base ring as the generator matrix.")

    size(v, 1) != 1 || return paritycheckmatrix(C) * v'
    return paritycheckmatrix(C) * v
end

function syndrome(v::Vector{Int64}, C::AbstractLinearCode)
    length(v) == length(C) ||
        error("Vector to be tested is of incorrect dimension; expected length $(C.n), received: $(size(v)).")
    return syndrome(matrix(field(C), v'), C)
end

"""
    in(v::Union{fq_nmod_mat, Vector{Int64}}, C::AbstractLinearCode)

Return whether or not `v` is a codeword of `C`.

The vector `v` is a valid codeword of `C` if and only if the syndrome of `v` is zero.
"""
in(v::Union{gfp_mat, fq_nmod_mat}, C::AbstractLinearCode) = iszero(syndrome(v, C))
in(v::Vector{Int64}, C::AbstractLinearCode) = iszero(syndrome(v, C))

"""
    ⊆(C1::AbstractLinearCode, C2::AbstractLinearCode)
    ⊂(C1::AbstractLinearCode, C2::AbstractLinearCode)
    issubcode(C1::AbstractLinearCode, C2::AbstractLinearCode)

Return whether or not `C1` is a subcode of `C2`.

A linear code `C1` is a subcode of another linear code `C2` if each row of
the generator matrix of `C1` is a valid codeword of `C2`.
"""
function ⊆(C1::AbstractLinearCode, C2::AbstractLinearCode)
    if field(C1) != field(C2) || length(C1) != length(C2) || dimension(C1) > dimension(C2)
        return false
    end

    # eachrow doesn't work on these objects
    for r in 1:dimension(C1)
        if C1.G[r, :] ∉ C2
            return false
        end
    end

    return true
end
⊂(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊆ C2
issubcode(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊆ C2

"""
    codecomplement(C1::AbstractLinearCode, C2::AbstractLinearCode)
    quo(C1::AbstractLinearCode, C2::AbstractLinearCode)
    quotient(C1::AbstractLinearCode, C2::AbstractLinearCode)
    /(C2::AbstractLinearCode, C1::AbstractLinearCode)

Return the linear code `C2 / C1` given `C1 ⊆ C2`.

Credit to Tommy Hofmann of the AbstractAlgebra/Nemo/Hecke packages for help with
debugging and providing the most elegant implementation now used here.
"""
function codecomplement(C1::AbstractLinearCode, C2::AbstractLinearCode)
    C1 ⊆ C2 || error("Need C1 ⊆ C2 in codecomplement(C1, C2) or C2 / C1.")
    F = field(C1)
    V = VectorSpace(F, length(C1))
    U, UtoV = sub(V, [V(C1.G[i, :]) for i in 1:nrows(C1.G)])
    W, WtoV = sub(V, [V(C2.G[i, :]) for i in 1:nrows(C2.G)])
    gensofUinW = [preimage(WtoV, UtoV(g)) for g in gens(U)]
    UinW, UinWtoW = sub(W, gensofUinW)
    Q, WtoQ = quo(W, UinW)
    C2modC1basis = [WtoV(x) for x in [preimage(WtoQ, g) for g in gens(Q)]]
    Fbasis = [[F(C2modC1basis[j][i]) for i in 1:dim(parent(C2modC1basis[1]))] for j in 1:length(C2modC1basis)]
    G = matrix(F, length(Fbasis), length(Fbasis[1]), vcat(Fbasis...))
    for r in 1:length(Fbasis)
        v = G[r, :]
        (v ∈ C2 && v ∉ C1) || error("Error in creation of basis for C2 / C1.")
    end
    return LinearCode(G)
end
quo(C1::AbstractLinearCode, C2::AbstractLinearCode) = codecomplement(C1, C2)
quotient(C1::AbstractLinearCode, C2::AbstractLinearCode) = codecomplement(C1, C2)
/(C2::AbstractLinearCode, C1::AbstractLinearCode) = codecomplement(C1, C2)

"""
    dual(C::AbstractLinearCode)

Return the dual of a linear code.

All of the necessary information for the dual is already stored in a LinearCode
object, so this implementation merely swaps data fields, e.g., `G <-> H`, without
doing any new computation.
"""
function dual(C::AbstractLinearCode)
    return LinearCode(field(C), length(C), length(C) - dimension(C), missing,
        deepcopy(paritycheckmatrix(C)), deepcopy(originalparitycheckmatrix(C)),
        deepcopy(generatormatrix(C)), deepcopy(originalgeneratormatrix(C)),
        deepcopy(paritycheckmatrix(C, true)), deepcopy(generatormatrix(C, true)),
        missing)
end

"""
    Hermitiandual(C::AbstractLinearCode)

Return the Hermitian dual of a linear code defined over a quadratic extension.
"""
Hermitiandual(C::AbstractLinearCode) = dual(LinearCode(Hermitianconjugatematrix(generatormatrix(C))))
# also valid to do LinearCode(Hermitianconjugatematrix(generatormatrix(dual(C))))

"""
    isequivalent(C1::AbstractLinearCode, C2::AbstractLinearCode)

Return whether or not `C1 ⊆ C2` and `C2 ⊆ C1`.
"""
isequivalent(C1::AbstractLinearCode, C2::AbstractLinearCode) = (C1 ⊆ C2 && C2 ⊆ C1)

"""
    isselfdual(C::AbstractLinearCode)

Return whether or not `isequivalent(C, dual(C))`.
"""
isselfdual(C::AbstractLinearCode) = isequivalent(C, dual(C))

"""
    isselforthogonal(C::AbstractLinearCode)
    isweaklyselfdual(C::AbstractLinearCode)

Return whether or not `C ⊆ dual(C)`.
"""
isselforthogonal(C::AbstractLinearCode) = C ⊆ dual(C)
isweaklyselfdual(C::AbstractLinearCode) = isselforthogonal(C)

"""
    ⊕(C1::AbstractLinearCode, C2::AbstractLinearCode)
    directsum(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊕ C2

Return the direct sum code of `C1` and `C2`.

The direct sum code has generator matrix `G1 ⊕ G2` and parity check matrix `H1 ⊕ H2`.
"""
function ⊕(C1::AbstractLinearCode, C2::AbstractLinearCode)
    field(C1) == field(C2) || error("Codes must be over the same field in the direct sum.")

    G = directsum(generatormatrix(C1), generatormatrix(C2))
    # this should just be the direct sum of these variables as well
    Gstand, Hstand = _standardform(G)
    return LinearCode(field(C1), length(C1), dimension(C1), missing, G, missing,
        directsum(C1.H, C2.H), missing, Gstand, Hstand, missing)
end
directsum(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊕ C2

"""
    ⊗(C1::AbstractLinearCode, C2::AbstractLinearCode)
    kron(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2
    tensorproduct(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2
    directproduct(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2
    productcode(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2

Return the (direct/tensor) product code of `C1` and `C2`.

The product code has generator matrix `G1 ⊗ G2` and parity check matrix `H1 ⊗ H2`.

# Notes
* The resulting product is not checked for any zero columns.
"""
function ⊗(C1::AbstractLinearCode, C2::AbstractLinearCode)
    field(C1) == field(C2) || error("Codes must be over the same field in the tensor product.")

    G = generatormatrix(C1) ⊗ generatormatrix(C2)
    Gstand, Hstand = _standardform(G)
    return LinearCode(field(C1), length(C1) * length(C2), dimension(C1) * dimension(C2),
        missing, G, missing, paritycheckmatrix(C1) ⊗ paritycheckmatrix(C2), missing,
        Gstand, Hstand, missing)
end
kron(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2
tensorproduct(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2
directproduct(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2
productcode(C1::AbstractLinearCode, C2::AbstractLinearCode) = C1 ⊗ C2

# # support is currently not implemented
# function support(C::AbstractLinearCode)
#     # returns the nonzero coefficient locations of the weight distribution
# end
#
# # support is currently not implemented
# function characteristicpolynomial(C::T) where T <: Union{LinearCode, CyclicCode}
#     # Return the characteristic polynomial of a linear code, as defined in [Lin1999]_.
#     R, x = PolynomialRing(Nemo.QQ, "x")
#     D = dual(C)
#     supD = support(D)
#     q = order(base_ring(C.G))
#     return q^(n - k) * prod([1 - x / j for j in supD if j > 0])
# end

"""
    extend(C::AbstractLinearCode)

Return the extended linear code of `C`.

This implementation chooses the most common form of extending a code, which is to
add an extra column to the generator matrix such that the sum of the coordinates
of each row is 0.
"""
function extend(C::AbstractLinearCode)
    col = [sum(C.G[i, :]) for i in 1:dimension(C)]
    M = MatrixSpace(field(C), dimension(C), 1)
    G = hcat(generatormatrix(C), M([characteristic(field(C)) - i for i in col]))
    M = MatrixSpace(field(C), 1, length(C) + 1)
    toprow = M([1 for _ in 1:(length(C) + 1)])
    M = MatrixSpace(field(C), length(C) - dimension(C), 1)
    rightcol = M(0)
    H = hcat(paritycheckmatrix(C), rightcol)
    H = vcat(toprow, H)
    Gstand, Hstand = _standardform(G)
    return LinearCode(field(C), size(G, 2), size(G, 1), missing,
        G, generatormatrix(C), H, paritycheckmatrix(C), Gstand, Hstand, missing)
end

"""
    puncture(C::AbstractLinearCode, cols::Vector{Int64})

Return the linear code of `C` punctured at the columns in `cols`.

Deletes the columns from the generator matrix and then removes any potentially
resulting zero rows.
"""
function puncture(C::AbstractLinearCode, cols::Vector{Int64})
    # to decode a punctured code ~kinda insert zeros into punctured locations
    # then uses original decoder + erasures

    !isempty(cols) || return C
    cols ⊆ 1:length(C) || error("Columns to puncture are not a subset of the index set.")
    length(cols) != length(C) || error("Cannot puncture all columns of a generator matrix.")

    G = C.G[:, setdiff(1:length(C), cols)]
    if rank(G) < size(G, 1)
        _, G = rref(G)
    end

    del = Vector{Int64}()
    for r in 1:size(G, 1)
        if iszero(G[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        G = G[setdiff(1:dimension(C), del), :]
    end
    _, H = right_kernel(G)
    # note the H here is transpose of the standard definition

    !(!iszero(G * H) || !iszero(H' * G')) ||
        error("Generator and parity check matrices are not transpose orthogonal.")

    Gstand, Hstand = _standardform(G)
    return LinearCode(field(C), size(G, 2), size(G, 1), missing, G, generatormatrix(C),
        H', paritycheckmatrix(C), Gstand, Hstand, missing)
end

"""
    expurgate(C::AbstractLinearCode, rows::Vector{Int64})

Return the linear code of `C` expuragated at the rows in `rows`.

Deletes the rows from the generator matrix and then removes any potentially
resulting zero columns.
"""
function expurgate(C::AbstractLinearCode, rows::Vector{Int64})
    !isempty(rows) || return C
    rows ⊆ 1:dimension(C) || error("Rows to expurgate are not a subset of the index set.")
    length(rows) != dimension(C) || error("Cannot expurgate all rows of a generator matrix.")

    G = C.G[setdiff(1:dimension(C), rows), :]
    # instead of calling LinearCode(G) here I want to keep C.G has Gorig so I am
    # repeating this constructor here such that I have more control over that
    if rank(G) < size(G, 1)
        _, G = rref(G)
    end

    del = Vector{Int64}()
    for r in 1:size(G, 1)
        if iszero(G[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        G = G[setdiff(1:dimension(C), del), :]
    end
    _, H = right_kernel(G)
    # note the H here is transpose of the standard definition

    !(!iszero(G * H) || !iszero(H' * G')) ||
        error("Generator and parity check matrices are not transpose orthogonal.")

    Gstand, Hstand = _standardform(G)
    return LinearCode(field(C), size(G, 2), size(G, 1), missing, G, generatormatrix(C),
        H', paritycheckmatrix(C), Gstand, Hstand, missing)
end

"""
    augment(C::AbstractLinearCode, M::fq_nmod_mat)

Return the linear code of `C` whose generator matrix is augmented with `M`.

Vertically joins the matrix `M` to the bottom of the generator matrix of `C`.
"""
function augment(C::AbstractLinearCode, M::Union{gfp_mat, fq_nmod_mat})
    length(C) == size(M, 2) || error("Rows to augment must have the same number of columns as the generator matrix.")
    field(C) == base_ring(M) || error("Rows to augment must have the same base field as the code.")

    G = vcat(generatormatrix(C), M)
    # instead of calling LinearCode(G) here I want to keep C.G has Gorig so I am
    # repeating this constructor here such that I have more control over that
    rk = rank(G)
    !iszero(rk) || error("Rank zero matrix passed into LinearCode constructor.")
    if rk < size(G, 1)
        _, G = rref(G)
    end

    del = Vector{Int64}()
    for r in 1:size(G, 1)
        if iszero(G[r, :])
            append!(del, r)
        end
    end

    # not necessary but faster to skip
    if !isempty(del)
        G = G[setdiff(1:dimension(C), del), :]
    end
    _, H = right_kernel(G)
    # note the H here is transpose of the standard definition

    !(!iszero(G * H) || !iszero(H' * G')) ||
        error("Generator and parity check matrices are not transpose orthogonal.")

    Gstand, Hstand = _standardform(G)
    for r in 1:size(Gstand, 1)
        iszero(Gstand[r, :] * H) || error("Column swap appeared in _standardform.")
    end

    return LinearCode(field(C), size(G, 2), size(G, 1), missing, G, generatormatrix(C),
        H', paritycheckmatrix(C), Gstand, Hstand, missing)
end

"""
    shorten(C::AbstractLinearCode, L::Vector{Int64})

Return the linear code of `C` shortened on the indices `L`.

Shortening is expurgating followed by puncturing. This implementation uses the
theorem that the dual of code shortened on `L` is equal to the puncture of the
dual code on `L`, i.e., `dual(puncture(dual(C), L))`.
"""
function shorten(C::AbstractLinearCode, L::Vector{Int64})
    # most parameter checks here done in puncture
    !isempty(L) || return C
    return dual(puncture(dual(C), L))
end

"""
    lengthen(C::AbstractLinearCode)

Return the lengthened linear code of `C`.

This augments the all 1's row and then extends.
"""
function lengthen(C::AbstractLinearCode)
    row = matrix(field(C), [1 for _ in 1:length(C)]')
    newcode = augment(C, row)
    return extend(newcode)
end

"""
    uuplusv(C1::AbstractLinearCode, C2::AbstractLinearCode, verify::Bool=true)
    Plotkinconstruction(C1::AbstractLinearCode, C2::AbstractLinearCode, verify::Bool=true)

Return the Plotkin- or so-called (u | u + v)-construction with `u ∈ C1` and `v ∈ C2`.

# Arguments
* `verify`: Runs a verification step on the sizes of the constructed code versus the
theoretically guarenteed values.
"""
function uuplusv(C1::AbstractLinearCode, C2::AbstractLinearCode, verify::Bool=true)
    # Returns the Plotkin `(u|u+v)`-construction with u = C1 and v = C2
    field(C1) == field(C2) || error("Base field must be the same in the Plotkin (u|u + v)-construction.")
    length(C1) == length(C2) || error("Both codes must have the same length in the Plotkin (u|u + v)-construction.")

    G = vcat(hcat(generatormatrix(C1), generatormatrix(C1)), hcat(parent(generatormatrix(C2))(0), generatormatrix(C2)))
    H = vcat(hcat(paritycheckmatrix(C1), parent(paritycheckmatrix(C1))(0)), hcat(-paritycheckmatrix(C2), paritycheckmatrix(C2)))
    Gstand, Hstand = _standardform(G)

    if ismissing(C1.d) || ismissing(C2.d)
        d = missing
    else
        d = minimum([2 * minimumdistance(C1), minimumdistance(C2)])
    end

    if verify
        size(G, 2) == 2 * length(C1) ||
            error("Something went wrong in the Plotkin (u|u + v)-construction;
                length ($(size(G, 2))) and expected length ($(2 * length(C1))) are not equal.")
        size(G, 1) == dimension(C1) + dimension(C2) ||
            error("Something went wrong in the Plotkin (u|u + v)-construction;
                dimension ($(size(G, 1))) and expected dimension ($(dimension(C1) + dimension(C2))) are not equal.")
    end

    return LinearCode(field(C1), size(G, 2), size(G, 1), d, G, missing, H, missing,
        Gstand, Hstand, missing)
end
Plotkinconstruction(C1::AbstractLinearCode, C2::AbstractLinearCode, verify::Bool=true) = uuplusv(C1, C2, verify)

"""
    subcode(C::AbstractLinearCode, k::Integer)

Return a `k`-dimension linear subcode code of `C`.
"""
function subcode(C::AbstractLinearCode, k::Integer)
    k >= 1 && k < dimension(C) || error("Cannot construct a $k-dimensional subcode of an $(dimension(C))-dimensional code.")
    k != dimension(C) || return C
    return LinearCode(generatormatrix(C)[1:k, :])
end

"""
    subcode(C::AbstractLinearCode, rows::Vector{Int64})

Return a linear subcode code of `C` using the rows of the generator matrix of `C`
listed in `rows`.
"""
function subcode(C::AbstractLinearCode, rows::Vector{Int64})
    !isempty(rows) || error("Row index set empty in subcode.")
    rows ⊆ 1:dimension(C) || error("Rows are not a subset of the index set.")
    length(rows) != dimension(C) || return C
    return LinearCode(generatormatri(C)[setdiff(1:dimension(C), rows), :])
end

"""
    juxtaposition(C1::AbstractLinearCode, C2::AbstractLinearCode)

Return the linear code generated by the horizontal concatenation of the generator
matrices of `C1` then `C2`.
"""
function juxtaposition(C1::AbstractLinearCode, C2::AbstractLinearCode)
    field(C1) == field(C2) || error("Cannot juxtapose two codes over different fields.")
    dimension(C1) == dimension(C2) || error("Cannot juxtapose two codes of different dimensions.")
    return LinearCode(hcat(generatormatrix(C1), generatormatrix(C2)))
end

"""
    constructionX(C1::AbstractLinearCode, C2::AbstractLinearCode, C3::AbstractLinearCode)

Return the linear code generated by the construction X procedure.

Let `C1` be an [n, k, d], `C2` be an [n, k - l, d + e], and `C3` be an [m, l, e] linear code
with `C2 ⊂ C1` be proper. Construction X creates a [n + m, k, d + e] linear code.
"""
function constructionX(C1::AbstractLinearCode, C2::AbstractLinearCode, C3::AbstractLinearCode)
    C1 ⊆ C2 || error("The first code must be a subcode of the second in construction X.")
    !isequivalent(C1, C2) || error("The first code must be a proper subcode of the second in construction X.")
    # the above line checks field(C1) == field(C2)
    field(C1) == field(C3) || error("All codes must be over the same base ring in construction X.")
    dimension(C2) == dimension(C1) + dimension(C3) ||
        error("The dimension of the second code must be the sum of the dimensions of the first and third codes.")
    # could do some verification steps on parameters
    return LinearCode(vcat(hcat(generatormatrix(C1 / C2), generatormatrix(C3)),
        hcat(generatormatrix(C1), zero_matrix(field(C1), dimension(C1), length(C3)))))
end

"""
    constructionX3(C1::AbstractLinearCode, C2::AbstractLinearCode, C3::AbstractLinearCode,
        C4::AbstractLinearCode, C5::AbstractLinearCode))

Return the linear code generated by the construction X3 procedure.

Let C1 = [n, k1, d1], C2 = [n, k2, d2], C3 = [n, k3, d3], C4 = [n4, k2 - k1, d4], and
C5 = [n5, k3 - k2, d5] with `C1 ⊂ C2 ⊂ C3`. Construction X3 creates an [n + n4 + n5, k3, d]
linear code with d ≥ min{d1, d2 + d4, d3 + d5}.
"""
function constructionX3(C1::AbstractLinearCode, C2::AbstractLinearCode, C3::AbstractLinearCode,
    C4::AbstractLinearCode, C5::AbstractLinearCode)

    C1 ⊆ C2 || error("The first code must be a subcode of the second in construction X3.")
    !isequivalent(C1, C2) || error("The first code must be a proper subcode of the second in construction X3.")
    C2 ⊆ C3 || error("The second code must be a subcode of the third in construction X3.")
    !isequivalent(C2, C3) || error("The second code must be a proper subcode of the third in construction X3.")
    # the above lines check field(C1) == field(C2) == field(C3)
    field(C1) == field(C4)  == field(C5) || error("All codes must be over the same base ring in construction X3.")
    dimension(C3) == dimension(C2) + dimension(C4) ||
        error("The dimension of the third code must be the sum of the dimensions of the second and fourth codes in construction X3.")
    dimension(C2) == dimension(C1) + dimension(C5) ||
        error("The dimension of the second code must be the sum of the dimensions of the first and fifth codes in construction X3.")

    C2modC1 = C2 / C1
    C3modC2 = C3 / C2
    F = field(C1)

    # could do some verification steps on parameters
    G = vcat(hcat(generatormatrix(C1), zero_matrix(F, dimension(C1), length(C4)), zero_matrix(F, dimension(C1), length(C5))),
    hcat(generatormatrix(C2modC1), generatormatrix(C4), zero_matrix(F, dimension(C4), length(C5))),
    hcat(generatormatrix(C3modC2), zero_matrix(F, dimension(C5), length(C4)), generatormatrix(C5)))
    return LinearCode(G)
end

"""
    upluswvpluswuplusvplusw(C1::AbstractLinearCode, C2::AbstractLinearCode)

Return the linear code generated by the (u + w | v + w | u + v + w)-construction.

Let C1 = [n, k1, d1] and C2 = [n, k2, d2]. This construction produces an [3n, 2k1 + k2]
linear code. For binary codes, wt(u + w | v + w | u + v + w) = 2 wt(u ⊻ v) - wt(w) + 4s,
where s = |{i | u_i = v_i = 0, w_i = 1}|.
"""
function upluswvpluswuplusvplusw(C1::AbstractLinearCode, C2::AbstractLinearCode)
    field(C1) == field(C2) || error("All codes must be over the same base ring in the (u + w | v + w | u + v + w)-construction.")
    length(C1) == length(C2) || error("All codes must be the same length in the (u + w | v + w | u + v + w)-construction.")

    # could do some verification steps on parameters
    G = vcat(hcat(generatormatrix(C1), zero_matrix(field(C1), dimension(C1), length(C1)), generatormatrix(C1)),
    hcat(generatormatrix(C2), generatormatrix(C2), generatormatrix(C2)),
    hcat(zero_matrix(field(C1), dimension(C1), length(C1)), generatormatrix(C1), generatormatrix(C1)))
    return LinearCode(G)
end

function expandedcode(C::AbstractLinearCode, K::FqNmodFiniteField, basis::Vector{fq_nmod})
    return LinearCode(expandmatrix(generatormatrix(C), K, basis))
end

# R.Pellikaan, On decoding by error location and dependent sets of error
# positions, Discrete Mathematics, 106107 (1992), 369-381.
# the Schur product of vector spaces is highly basis dependent and is often the
# full vector space (an [n, n, 1] code)
function entrywiseproductcode(C::AbstractLinearCode, D::AbstractLinearCode)
    field(C) == field(D) || error("Codes must be over the same field in the Schur product.")
    length(C) == length(D) || error("Codes must have the same length in the Schur product.")

    GC = generatormatrix(C)
    GD = generatormatrix(D)
    indices = Vector{Tuple{Int, Int}}()
    for i in 1:dimension(C)
        for j in 1:dimension(D)
            i <= j && push!(indices, (i, j))
        end
    end
    return LinearCode(matrix(field(C), vcat([GC[i, :] .* GD[j, :] for (i, j) in indices]...)))

    # verify C ⊂ it
end
*(C::AbstractLinearCode, D::AbstractLinearCode) = entrywiseproductcode(C, D)
Schurproductcode(C::AbstractLinearCode, D::AbstractLinearCode) = entrywiseproductcode(C, D)
Hadamardproductcode(C::AbstractLinearCode, D::AbstractLinearCode) = entrywiseproductcode(C, D)
componentwiseproductcode(C::AbstractLinearCode, D::AbstractLinearCode) = entrywiseproductcode(C, D)

function VectorSpace(C::AbstractLinearCode)
    V = VectorSpace(field(C), length(C))
    G = generatormatrix(C)
    return sub(V, [V(G[i, :]) for i in 1:nrows(G)])
end

# needs significant testing, works in a very special case
function evensubcode(C::AbstractLinearCode)
    F = field(C)
    # add checks here that this is only implemented for binary codes

    VC, ψ = VectorSpace(C)
    GFVS = VectorSpace(F, 1)
    homo1quad = ModuleHomomorphism(VC, GFVS, matrix(F, dim(VC), 1, vcat([F(weight(ψ(g).v) % 2) for g in gens(VC)]...)))
    evensub, ϕ1 = kernel(homo1quad)
    !iszero(dim(evensub)) && return LinearCode(vcat([ψ(ϕ1(g)).v for g in gens(evensub)]...))
    return missing
end

# needs significant testing, works in a very special case
function doublyevensubcode(C::AbstractLinearCode)
    F = field(C)
    # add checks here that this is only implemented for binary codes

    VC, ψ = VectorSpace(C)
    GFVS = VectorSpace(F, 1)

    # first get the even subspace
    homo1quad = ModuleHomomorphism(VC, GFVS, matrix(F, dim(VC), 1, vcat([F(weight(ψ(g).v) % 2) for g in gens(VC)]...)))
    evensub, ϕ1 = kernel(homo1quad)

    if !iszero(dim(evensub))
        # now control the overlap (Ward's divisibility theorem)
        homo2bi = ModuleHomomorphism(evensub, evensub, matrix(F, dim(evensub), dim(evensub),
            vcat([F(weight(matrix(F, 1, length(C), ψ(ϕ1(gens(evensub)[i])).v .* ψ(ϕ1(gens(evensub)[j])).v)) % 2)
            for i in 1:dim(evensub), j in 1:dim(evensub)]...)))
        evensubwoverlap, μ1 = kernel(homo2bi)

        if !iszero(dim(evensubwoverlap))
            # now apply the weight four condition
            homo2quad = ModuleHomomorphism(evensubwoverlap, GFVS, matrix(F, dim(evensubwoverlap), 1,
                vcat([F(div(weight(ψ(ϕ1(μ1(g))).v), 2) % 2) for g in gens(evensubwoverlap)]...)))
            foursub, ϕ2 = kernel(homo2quad)

            if !iszero(dim(foursub))
                return LinearCode(vcat([ψ(ϕ1(μ1(ϕ2(g)))).v for g in gens(foursub)]...))
            end
        end
    end
    return missing
end

# currently incomplete
# function triplyevensubcode(C::AbstractLinearCode)
#     F = field(C)
#     VC, ψ = VectorSpace(C)
#     GFVS = VectorSpace(F, 1)
#
#     # first get the even subspace
#     homo1quad = ModuleHomomorphism(VC, GFVS, matrix(F, dim(VC), 1, vcat([F(weight(ψ(g).v) % 2) for g in gens(VC)]...)))
#     evensub, ϕ1 = kernel(homo1quad)
#
#     # now control the overlap (Ward's divisibility theorem)
#     homo2bi = ModuleHomomorphism(evensub, evensub, matrix(F, dim(evensub), dim(evensub),
#         vcat([F(weight(matrix(F, 1, length(C), ψ(ϕ1(gens(evensub)[i])).v .* ψ(ϕ1(gens(evensub)[j])).v)) % 2)
#         for i in 1:dim(evensub), j in 1:dim(evensub)]...)))
#     evensubwoverlap, μ1 = kernel(homo2bi)
#
#     # now apply the weight four condition
#     homo2quad = ModuleHomomorphism(evensubwoverlap, GFVS, matrix(F, dim(evensubwoverlap), 1,
#         vcat([F(div(weight(ψ(ϕ1(μ1(g))).v), 2) % 2) for g in gens(evensubwoverlap)]...)))
#     foursub, ϕ2 = kernel(homo2quad)
#
#     # now control the triple overlap
#     ###########
#     # how do we actually represent this?
#
#
#     # now apply the weight eight condition
#     homo3 = ModuleHomomorphism(foursub, GFVS, matrix(F, dim(foursub), 1, vcat([F(div(weight(ψ(ϕ1(ϕ2(g))).v), 4) % 2) for g in gens(foursub)]...)))
#     eightsub, ϕ3 = kernel(homo3)
#     # # println(dim(eightsub))
#     # # println(vcat([ψ(ϕ1(ϕ2(ϕ3(g)))).v for g in gens(eightsub)]...))
#     #
#     # q, γ = quo(foursub, eightsub)
#     # println(dim(q))
#     # # println([preimage(γ, g).v for g in gens(q)])
#     # return LinearCode(vcat([ψ(ϕ1(ϕ2(preimage(γ, g)))).v for g in gens(q)]...))
#     # # return LinearCode(vcat([ψ(ϕ1(ϕ2(ϕ3(g)))).v for g in gens(eightsub)]...))
#
#     return LinearCode(vcat([ψ(ϕ1(μ1(ϕ2(g)))).v for g in gens(foursub)]...))
# end

max_stars_repo_name: esabo/CodingTheory

max_issues_repo_name: esabo/CodingTheory

max_forks_repo_name: esabo/CodingTheory


text:
for op in (:check, :save, :load, :delete)
    eval(quote
        function $op(d::DataEntry, v, file_name::String, groups...)
            folder_path = get_folder_path(d, v)
            file_path = joinpath(folder_path, file_name)
            ext = _get_extension(file_name)
            if isempty(groups) || ext == "csv"
                group_path = ""
            else
                group_path = get_group_path(d, v)
            end
            $op(file_path, group_path, ext, groups...)
        end
    end)
end

function check(file_path::AbstractString, group_path::AbstractString, ext, groups...)
    if isfile(file_path)
        f = _get_function("check_", ext)
        all([f(file_path, group_path, group_name) for group_name in groups])
    else
        false
    end
end

"""
$(SIGNATURES)

Save data in `groups` in to file `file_path`.
"""
function save(file_path::AbstractString, group_path::AbstractString, ext, groups...)
    f = _get_function("save_", ext)
    f(file_path, group_path, groups...)
end

function load(file_path::AbstractString, group_path::AbstractString, ext, groups...)
    _check_file(file_path)
    f = _get_function("load_", ext)
    f(file_path, group_path, groups...)
end

function delete(
    file_path::AbstractString,
    group_path::AbstractString,
    ext,
    groups...,
)
    if isempty(groups)
        rm(file_path)
    elseif isfile(file_path)
        f = _get_function("delete_", ext)
        f(file_path, group_path, groups...)
    end
end

function load_file_array(d::DataEntry, v, file_name::String, groups...)
    r_str = Regex(file_name * "-(\\d+)\\.(.+)")
    folder_path = get_folder_path(d, v)
    f_list = readdir(folder_path)
    group_path = get_group_path(d, v)
    # names = [group_path * "/" * i for i in groups]
    res = []
    for f in f_list
        m = match(r_str, f)
        if m != nothing
            f_path = joinpath(folder_path, m.match)
            ext = _get_extension(m.match)
            push!(res, load(f_path, group_path, ext, groups...))
        end
    end
    res
end

ext: jl


text:
"""
Run tests from fmi-cross-check Git repository.
"""

using Test
using LibGit2

thisDir = dirname(Base.source_path())
fmiCrossCheckDir = joinpath(thisDir,"fmi-cross-check")

include("$(dirname(thisDir))/src/FMUSimulator.jl")

if Sys.iswindows()
    fmiCrossCheckFMUDir = joinpath(thisDir, "fmi-cross-check", "fmus", "2.0", "me", "win$(Sys.WORD_SIZE)")
elseif Sys.islinux()
    fmiCrossCheckFMUDir = joinpath(thisDir, "fmi-cross-check", "fmus", "2.0", "me", "linux$(Sys.WORD_SIZE)")
else
    error("OS not supportet for this tests.")
end

"""
    updateFMICrossTest()

Clone or fetch and pull modelica/fmi-cross-check repository from
https://github.com/modelica/fmi-cross-check.git.
"""
function updateFmiCrossCheck()

    if isdir(fmiCrossCheckDir)
        # Update repository
        println("Updating repository modelica/fmi-cross-check.")
        repo = GitRepo(fmiCrossCheckDir)
        LibGit2.fetch(repo)
        LibGit2.merge!(repo, fastforward=true)

    else
        # Clone repository
        println("Cloning repository modelica/fmi-cross-check.")
        println("This may take some minutes.")
        LibGit2.clone("https://github.com/modelica/fmi-cross-check.git", fmiCrossCheckDir )
        println("Cloned modelica/fmi-cross-check successfully.")
    end
end

"""
    cleanFmiCrossCheck()

Resets fmi-cross-check repository. All changes will get lost!
"""
function cleanFmiCrossCheck()

    if isdir(fmiCrossCheckDir)
        repo = GitRepo(fmiCrossCheckDir)
        head_oid = LibGit2.head_oid(repo)
        mode = LibGit2.Consts.RESET_HARD
        LibGit2.reset!(repo, head_oid, mode)
    end
end

"""
    runFMICrossTests()

Run and test all fmi-cross-test that are supportet on current system.
"""
function runFMICrossTests()

    # Check if repository is up to date
    updateFmiCrossCheck()

    @testset "FMI Cross Check" begin
        # Windows tests
        if Sys.iswindows()
            @testset "Catia" begin
                testCatiaFMUs()
            end;
            @testset "SystemModeler" begin
                testSystemModelerFMUs()
            end;
            @testset "Test-FMUs" begin
                testTestFMUsFMUs()
            end;

        # Linux tests
        elseif Sys.islinux()
            @testset "JModelica.org" begin
                testJModelicaFMUs()
            end;
            @testset "MapleSim" begin
                testMapleSimFMUs()
            end;
            @testset "SystemModeler" begin
                testSystemModelerFMUs()
            end;
            @testset "SystemModeler" begin
                testSystemModelerFMUs()
            end;
            @testset "Test-FMUs" begin
                testTestFMUsFMUs()
            end;
        end
    end;
end

"""
    testCatiaFMUs()

Test ModelExchange 2.0 FMUs generated by Catia

Tests versions R2015x and R2016x on 32 and 64 bit Windows.

# Tests
* modelBooleanNetwork1
* ControlledTemperature
* CoupledClutches
* DFFREG
* MixtureGases
* Rectifier
"""
function testCatiaFMUs()

    if !Sys.iswindows()
        error("Test only supportet on 32 and 64 bit Windows.")
    end

    toolName = "CATIA"
    versions = ["R2015x", "R2016x"]
    tests = ["BooleanNetwork1" "ControlledTemperature" "CoupledClutches" "DFFREG" "" "Rectifier";
             "BooleanNetwork1" "ControlledTemperature" "CoupledClutches" "DFFREG" "MixtureGases" "Rectifier"]

    testTool(toolName, versions, tests)
end

"""
    testJModelicaFMUs()

Test ModelExchange 2.0 FMUs generated by JModelica.org.

Tests version 1.15 on 64 bit Linux.

# Tests
* CauerLowPassAnalog
* ControlledTemperature
* CoupledClutches
* PID_Controller
"""
function testJModelicaFMUs()

    if !Sys.islinux() || Sys.WORD_SIZE != 64
        error("Test only supportet on 64 bit Linux.")
    end

    toolName = "JModelica.org"
    versions = ["1.15"]
    tests = ["CauerLowPassAnalog" "ControlledTemperature" "CoupledClutches" "PID_Controller"]

    testTool(toolName, versions, tests)
end

"""
    testMapleSimFMUs()

Test ModelExchange 2.0 FMUs generated by MapleSim.

Tests versions 7.01, 2015.1, 2015.2, 2016.1, 2016.2, 2018 on 32 or 64 bit Linux.

# Tests
* ControlledTemperature
* CoupledClutches
* Rectifier
"""
function testMapleSimFMUs()

    if !Sys.islinux()
        error("Test only supportet on 32 and 64 bit Linux.")
    end

    toolName = "MapleSim"
    if Sys.WORD_SIZE == 32
         versions = ["7.01", "2015.1", "2015.2"]
     elseif Sys.WORD_SIZE == 64
         versions = ["7.01", "2015.1", "2015.2", "2016.1", "2016.2", "2018"]
    else
        error("Unknown WORD_SIZE: $WORD_SIZE")
     end

    tests = ["ControlledTemperature" "" "";
             "ControlledTemperature" "CoupledClutches" "Rectifier";
             "ControlledTemperature" "CoupledClutches" "Rectifier";
             "ControlledTemperature" "CoupledClutches" "Rectifier";
             "ControlledTemperature" "CoupledClutches" "Rectifier";
             "ControlledTemperature" "CoupledClutches" "Rectifier"]

    testTool(toolName, versions, tests)
end

"""
    testSystemModelerFMUs()

Test ModelExchange 2.0 FMUs generated by SystemModeler.

Tests version 5.0 on 64 bit Linux or Windows.

# Tests
* ControlledTemperature
* CoupledClutches
"""
function testSystemModelerFMUs()

    if Sys.WORD_SIZE != 64
        error("Test only supportet on 64 bit Linux or Windows.")
    end

    toolName = "SystemModeler"
    versions = ["5.0"]

    tests = ["ControlledTemperature" "CoupledClutches"]

    testTool(toolName, versions, tests)
end

"""
    testTestFMUsFMUs()

Test ModelExchange 2.0 FMUs generated by Test-FMUs.

Tested version 0.0.1 on 64 bit Linux or Windows

# Tests
*
"""
function testTestFMUsFMUs()

    if Sys.WORD_SIZE != 64
        error("Test only supportet on 64 bit Linux or Windows.")
    end

    toolName = "Test-FMUs"
    versions = ["0.0.1"]

    tests = ["BouncingBall" "Dahlquist" "Feedthrough" "Resource" "Stair" "VanDerPol"]

    testTool(toolName, versions, tests)
end

"""
    function testTool(toolName::String, versions::Array{String,1}, tests)

Heler function to test for generic tools, versions and test cases.
"""
function testTool(toolName::String, versions::Array{String,1}, tests)

    for (i,version) in enumerate(versions)
        @testset "$version" begin
            for test in tests[i,:]
                if test != ""
                    @test begin
                        model = joinpath(fmiCrossCheckFMUDir, "$toolName", "$version", "$test", "$test.fmu")
                        main(model)
                    end;
                end
            end
        end;
    end
end

max_stars_repo_name: cuihantao/JuliaFMI

max_issues_repo_name: cuihantao/JuliaFMI

max_forks_repo_name: cuihantao/JuliaFMI


text:
@test rowcoord(ZIPCODES[1, :]) == (43.005895,-71.013202)
@test length(rowcoord(ZIPCODES)) == nrow(ZIPCODES)


text:
# funcs.jl - common calculation functions

export interest, rate, mortgage

function interest(amount, rate)
    return amount * (1 + rate)
end

function rate(amount, interest)
    return interest / amount
end

# uses payment function from Mortgage.jl
function mortgage(home_price, down_payment, rate, years)
    return payment(home_price - down_payment, rate, years)
end


text:
@testset "latex printer" begin

    tests = [
        ("a(1)", "a_{t+1}",),
        ("astar(1)", "a_{t+1}^{\\star}",),
        ("a_i_j__k(1)", "a_{i,j,t+1}^{k}",),
        ("a_i_j__k(+1)", "a_{i,j,t+1}^{k}",),
        ("a_i_j__k(-1)", "a_{i,j,t-1}^{k}",),
        ("a__k_i_j(1)", "a_{i,j,t+1}^{k}",),
        ("abar_xbar(1)", "\\overline{a}_{\\overline{x},t+1}"),
        ("adot_xbar(1)", "\\dot{a}_{\\overline{x},t+1}"),
        ("abar_xdot(1)", "\\overline{a}_{\\dot{x},t+1}"),
        ("adot_xdot(1)", "\\dot{a}_{\\dot{x},t+1}"),
        ("abar_i_jstar__k(1)", "\\overline{a}_{i,j^{\\star},t+1}^{k}"),
        ("a__p1_p2(1)", "a_{p2,t+1}^{p1}"),
        ("(1 - delta) * k", "\\left(1-\\delta\\right) k")
    ]

    for (have, want) in tests
        @test Dolang.latex(have) == want
    end

end  # testset

max_issues_repo_name: EconForge/Dolang


text:
function KS_solve_Emin_cg( LF::LF3dGrid, Gv::GvectorsT,
                          V_ionic, Focc, Nstates::Int64;
                          v0 = nothing,
                          Potentials0 = nothing,
                          α_t=3e-5, Niter=1000,
                          verbose=false )
    #
    Npoints = LF.Nx * LF.Ny * LF.Nz
    ΔV = LF.LFx.h * LF.LFy.h * LF.LFz.h
    #
    # Setup trial solution
    #
    if v0 == nothing
        # random trial solution
        srand(1234)
        v = rand( Npoints, Nstates )
        v = orthonormalize( LF, v )
    else
        v = copy(v0)
        # v0 is assumed to be ortonormalized
    end
    #
    # Setup potentials
    #
    if Potentials0 == nothing
        # initialize potential: Ionic, Hartree, XC
        Potentials = PotentialsT( V_ionic, zeros(Npoints), zeros(Npoints) )
    else
        Potentials = copy(Potentials0)
    end

    # Initialize energy components: Total, Kinetic, Ionic, Hartree, XC
    Energies = EnergiesT( 0.0, 0.0, 0.0, 0.0, 0.0 )
    #
    # Calculate energies and update potetials
    #
    calc_Energies!( LF, Gv, Energies, Potentials, Focc, v )
    #
    Etot_old = Energies.Total
    print_Energies( Energies )
    #
    g     = zeros( Float64, Npoints, Nstates )
    g_old = zeros( Float64, Npoints, Nstates )
    d     = zeros( Float64, Npoints, Nstates )
    d_old = zeros( Float64, Npoints, Nstates )
    #
    β = 0
    α = 0
    for iter = 1:Niter
        #
        g = calc_grad( LF, Potentials, Focc, v )
        #
        if iter != 1
            β = trace( g' * g)/trace( g_old'*g_old )
            #β = trace( (g-g_old)'*g ) / trace( g_old'*g_old )
            #β = trace( (g-g_old)'*g )/ trace( (g-g_old)'*d )
        end
        d = -g + β*d_old
        #
        # compute gradient at trial step
        #
        v2 = orthonormalize(LF, v + α_t*d)
        rho = calc_rho( Focc, v2 )
        V_Hartree = solve_poisson_FFT( Gv, rho )
        V_xc = excVWN( rho ) + rho .* excpVWN( rho )
        wPot = PotentialsT( Potentials.Ionic, V_Hartree, V_xc )
        g_t = calc_grad( LF, wPot, Focc, v2 )
        #
        # compute estimate of best step and update current trial vectors
        #
        denum = trace( (g - g_t)' * d )
        if denum != 0.0
            α = abs( α_t * trace( g'*d )/denum )
        else
            α = 0.0
        end
        #
        v = v + α*d
        #
        # Calculate energies and update potetials
        #
        v = orthonormalize( LF, v )
        calc_Energies!( LF, Gv, Energies, Potentials, Focc, v )
        Etot = Energies.Total
        #
        if verbose
            @printf("%8d %16.9f %16.9e\n", iter, Etot, abs(Etot-Etot_old))
        end
        if abs(Etot-Etot_old) < 1.e-7
            if verbose
                @printf("Emin CG converges in %8d iterations\n", iter)
            end
            break
        end
        Etot_old = Etot
        #
        g_old = g[:,:]
        d_old = d[:,:]
    end
    return Energies, v, Potentials
end


text:
function solve(
    i::BipartiteMatchingInstance{T},
    ::DynamicProgramming;
    kwargs...,
) where {T}
    return matching_dp(i; kwargs...)
end

function matching_dp_imperfect(instance::BipartiteMatchingInstance{T}) where {T}
    state = Dict{Tuple{Int, Int}, Float64}(
        (i, j) => -Inf for i in instance.vertex_left,
        j in instance.vertex_right
    )
    solutions = Dict{Tuple{Int, Int}, Vector{Edge{T}}}()

    # Initialise: j=1 and i=1.
    v1 = first(instance.vertex_right)
    for i in instance.vertex_left
        if has_edge(instance.graph, i, v1)
            state[i, v1] = instance.reward[Edge(i, v1)]
            solutions[i, v1] = [Edge(i, v1)]
        end
    end

    v1 = first(instance.vertex_left)
    for j in instance.vertex_right
        if has_edge(instance.graph, v1, j)
            state[v1, j] = instance.reward[Edge(v1, j)]
            solutions[v1, j] = [Edge(v1, j)]
        end
    end

    # Main iteration, considering the first vertex of each side is already done.
    for (i, prev_i) in
        zip(instance.vertex_left[2:end], instance.vertex_left[1:(end - 1)])
        for (j, prev_j) in zip(
            instance.vertex_right[2:end],
            instance.vertex_right[1:(end - 1)],
        )
            # Take (i, j).
            if has_edge(instance.graph, i, j) && state[prev_i, prev_j] > -Inf
                state[i, j] =
                    instance.reward[Edge(i, j)] + state[prev_i, prev_j]
                solutions[i, j] = solutions[prev_i, prev_j]
                push!(solutions[i, j], Edge(i, j))
            end

            # Don't take (i, j): copy a previous solution.
            if state[prev_i, j] > state[i, j]
                state[i, j] = state[prev_i, j]
                solutions[i, j] = solutions[prev_i, j]
            end
            if state[i, prev_j] > state[i, j]
                state[i, j] = state[i, prev_j]
                solutions[i, j] = solutions[i, prev_j]
            end
        end
    end

    val = state[last(instance.vertex_left), last(instance.vertex_right)]
    sol = solutions[last(instance.vertex_left), last(instance.vertex_right)]

    return BipartiteMatchingSolution(instance, sol, val)
    # return state, solutions
end

function matching_dp(instance::BipartiteMatchingInstance{T}) where {T}
    # First, determine the size of a perfect matching.
    nr = Dict(k => 1.0 for (k, v) in instance.reward)
    ni = BipartiteMatchingInstance(instance, nr)
    ns = matching_dp_imperfect(ni)
    m = length(ns.solution)

    # If solving the matching as imperfect gives the right size, then done.
    naive_sol = matching_dp_imperfect(instance)
    if length(naive_sol.solution) == m
        return naive_sol
    end

    # Otherwise, solve this matching to ensure the perfect constraint is met.
    # TODO: reuse the state from the previous step?

    # Three indices: two vertices, like matching_dp_imperfect; the number of edges still to add to have a perfect matching.
    state = Dict{Tuple{Int, Int, Int}, Float64}(
        (i, j, k) => -Inf for i in instance.vertex_left,
        j in instance.vertex_right, k in 0:m
    )
    solutions = Dict{Tuple{Int, Int, Int}, Vector{Edge{T}}}(
        (i, j, m) => Edge{T}[] for i in instance.vertex_left,
        j in instance.vertex_right
    )

    # Initialise: j=1 and i=1, k=m-1.
    v1 = first(instance.vertex_right)
    for i in instance.vertex_left
        if has_edge(instance.graph, i, v1)
            state[i, v1, m - 1] = instance.reward[Edge(i, v1)]
            solutions[i, v1, m - 1] = [Edge(i, v1)]
        end
    end

    v1 = first(instance.vertex_left)
    for j in instance.vertex_right
        if has_edge(instance.graph, v1, j)
            state[v1, j, m - 1] = instance.reward[Edge(v1, j)]
            solutions[v1, j, m - 1] = [Edge(v1, j)]
        end
    end

    # Main iteration, considering the first vertex of each side is already done.
    for k in (m - 2):-1:0
        for (i, prev_i) in
            zip(instance.vertex_left[2:end], instance.vertex_left[1:(end - 1)])
            for (j, prev_j) in zip(
                instance.vertex_right[2:end],
                instance.vertex_right[1:(end - 1)],
            )
                # Take (i, j).
                if has_edge(instance.graph, i, j) &&
                   state[prev_i, prev_j, k + 1] > -Inf
                    state[i, j, k] =
                        instance.reward[Edge(i, j)] +
                        state[prev_i, prev_j, k + 1]
                    solutions[i, j, k] = solutions[prev_i, prev_j, k + 1]
                    push!(solutions[i, j, k], Edge(i, j))
                end

                # Don't take (i, j): mix previous solutions (as a special case, maybe take a previous solution completely).
                for k2 in 0:m
                    for k3 in 0:m
                        if state[prev_i, j, k2] + state[i, prev_j, k3] >
                           state[i, j, k]
                            potential_sol = union(
                                Set{Edge{T}}(solutions[prev_i, j, k2]),
                                Set{Edge{T}}(solutions[i, prev_j, k3]),
                            )
                            val = sum(
                                instance.reward[Edge(src(e), dst(e))] for
                                e in potential_sol
                            )

                            if length(potential_sol) == m - k &&
                               val > state[i, j, k]
                                state[i, j, k] = val
                                solutions[i, j, k] = collect(potential_sol)
                            end
                        end
                    end
                end
            end
        end
    end

    val = state[last(instance.vertex_left), last(instance.vertex_right), 0]
    sol = solutions[last(instance.vertex_left), last(instance.vertex_right), 0]

    return BipartiteMatchingSolution(instance, sol, val)
end


text:
#Artificial viscosity

#Artificial viscosity according to Colella & Woodward 1984 (eq. 4.5)
function avisc_CW(u, v, g, nx, ny, dx, dy)

    const cvisc = 0.1

    avisco_x = zeros(ny, nx)
    avisco_y = zeros(ny, nx)

    for j = (g-1):(ny-g+1), i = (g-1):(nx-g+1)

        #x-interface
        divU_x = (u[j, i] - u[j, i-1])/dx + 0.25(v[j+1, i] + v[j+1, i-1] - v[j-1, i] - v[j-1, i-1])/dx
        avisco_x[j, i] = cvisc*max(-divU_x*dx, 0.0)

        #y-interface value
        divU_y = 0.25(u[j, i+1] + u[j-1, i+1] - u[j, i-1] - u[j-1, i-1])/dx + (v[j, i] - v[j-1, i])/dy
        avisco_y[j, i] = cvisc*max(-divU_y*dy, 0.0)
    end

    return avisco_x, avisco_y
end

#General function for artificial viscosity
function artificial_viscosity(hyd::data2d, dt)

    dx = abs(hyd.x[2]-hyd.x[1])
    dy = abs(hyd.y[2]-hyd.y[1])
    aviscx, aviscy = avisc_CW(hyd.velx, hyd.vely, hyd.g, hyd.nx, hyd.ny, dx, dy)

    fx = zeros(hyd.ny, hyd.nx, 4)
    fy = zeros(hyd.ny, hyd.nx, 4)

    for j = (hyd.g):(hyd.ny-hyd.g+1), i = (hyd.g):(hyd.nx-hyd.g+1)
        for k = 1:4
            fx[j, i, k] = aviscx[j, i]*(hyd.q[j, i-1, k] - hyd.q[j, i, k])
            fy[j, i, k] = aviscy[j, i]*(hyd.q[j-1, i, k] - hyd.q[j, i, k])
        end
    end

    return fx, fy
end

max_stars_repo_name: natj/hydro

max_issues_repo_name: natj/hydro

max_forks_repo_name: natj/hydro


text:
export DataDrill, DataDrillStartOnly, zchars, DataDrillPrimitive

# Types to define Initial vs Development Drilling
#------------------------------------------

abstract type AbstractRegimeType end
struct InitialDrilling     <: AbstractRegimeType end
struct DevelopmentDrilling <: AbstractRegimeType end
struct FinishedDrilling    <: AbstractRegimeType end

+(::InitialDrilling, i) = DevelopmentDrilling()
+(::DevelopmentDrilling, i) = FinishedDrilling()
+(::FinishedDrilling, i) = nothing

-(::InitialDrilling, i) = nothing
-(::DevelopmentDrilling, i) = InitialDrilling()
-(::FinishedDrilling, i) = DevelopmentDrilling()

==(::AbstractRegimeType, i) = false
==(::A, ::A) where {A<:AbstractRegimeType} = true

isless(::A, ::A) where {A<:AbstractRegimeType} = false
isless(::InitialDrilling,  ::Union{DevelopmentDrilling,FinishedDrilling}) = true
isless(::FinishedDrilling, ::Union{InitialDrilling,DevelopmentDrilling})  = false
isless(::DevelopmentDrilling, ::InitialDrilling)  = false
isless(::DevelopmentDrilling, ::FinishedDrilling) = true

# DataSet
#---------------------------

function x_in_statespace(x, wp::AbstractUnitProblem)
    setx = Set(x)
    nS = length(wp)
    if all(x .== 0)
        check = true
    else
        check = issubset( setx, 1:nS )
    end
    check || throw(error("not all x ∈ $setx are in $wp"))
    return check
end

x_in_statespace(x, wp) = true

function DataDrillChecks(j1ptr, j2ptr, tptr, jtstart, jchars, ichars, y, x, zchars, wp)

    # check i
    length(j1ptr)-1 == length(ichars) == length(j2ptr)   ||
        throw(DimensionMismatch("Lengths of ichars, j2ptr, and j1ptr must agree"))

    # check j
    last(j1ptr) == first(j2ptr)                          ||
        throw(error("j1ptr and j2ptr must be conts"))

    # check t
    last(j2ptr) == length(jtstart) == length(tptr)-1     ||
        throw(DimensionMismatch("lengths of tptr, jtstart must be consistent"))

    # check tchars
    if length(y) == length(x) == last(tptr)-1
    else
        leny = length(y)
        lenx = length(x)
        tlast = last(tptr)-1
        throw(DimensionMismatch("length(y) = $leny, length(x) = $lenx, last(tptr)-1=$tlast"))
    end
    # pointers are sorted
    issorted(j1ptr) && issorted(j2ptr) && issorted(tptr) ||
        throw(error("Pointers not sorted"))

    # time vars are OK
    for j in 1:length(tptr)-1
        0 < jtstart[j] || throw("jtstart[$j] = $(jtstart[j]) <= 0")
        jtstart[j] + tptr[j+1] - 1 - tptr[j] <= length(zchars) ||
            throw(error("don't have z for all times implied by jtstart"))
    end

    x_in_statespace(x,wp)

    return true
end

struct DataDrillPrimitive{
        R<:AbstractStaticPayoff, ETV<:ExogTimeVars,
        ITup<:Tuple, XT, UP<:AbstractStateSpace
        } <: AbstractDataDrill
    reward::R

    j1ptr::Vector{Int}    # tptr, jchars is in  j1ptr[i] : j1ptr[i+1]-1
    j2ptr::UnitRange{Int} # tptr, jchars is in  j2ptr[i]
    tptr::Vector{Int}     # tchars       is in  tptr[j] : tptr[j+1]-1
    jtstart::Vector{Int}  # zvars for lease j start at zchars[jtstart(data,j)]

    # leases per unit
    j1chars::Vector{Float64} # weights for lease observations

    # drilling histories
    ichars::Vector{ITup}
    y::Vector{Int}
    x::Vector{XT}

    # time indices
    zchars::ETV
    wp::UP

    function DataDrillPrimitive(
        reward::R, j1ptr, j2ptr, tptr, jtstart,
        jchars, ichars::Vector{ITup}, y, x::Vector{XT}, zchars::ETV, wp::UP
    ) where {
        R, ETV, ITup, XT, UP
    }

        chk = DataDrillChecks(j1ptr, j2ptr, tptr, jtstart, jchars, ichars, y, x, zchars, wp)
        if chk == false
            throw(error("data didn't check out!"))
        end

        return new{R,ETV,ITup,XT,UP}(
            reward, j1ptr, j2ptr, tptr, jtstart, jchars, ichars, y, x, zchars, wp
        )
    end
end

reward(x::DataDrillPrimitive) = x.reward

struct DataDrill{M<:AbstractDrillModel, ETV<:ExogTimeVars, ITup<:Tuple, XT} <: AbstractDataDrill
    model::M

    j1ptr::Vector{Int}             # tptr, jchars is in  j1ptr[i] : j1ptr[i+1]-1
    j2ptr::UnitRange{Int}          # tptr, jchars is in  j2ptr[i]
    tptr::Vector{Int}              # tchars       is in  tptr[j] : tptr[j+1]-1
    jtstart::Vector{Int}           # zvars for lease j start at zchars[jtstart(data,j)]

    # leases per unit
    j1chars::Vector{Float64}       # weights for lease observations

    # drilling histories
    ichars::Vector{ITup}
    y::Vector{Int}
    x::Vector{XT}

    # time indices
    zchars::ETV
    dtv::DrillingTmpVars{Float64}

    function DataDrill(model::M, j1ptr, j2ptr, tptr, jtstart,
        jchars, ichars::Vector{ITup}, y, x::Vector{XT}, zchars::ETV
    ) where {M, ETV, ITup, XT}

        wp = statespace(model)
        chk = DataDrillChecks(j1ptr, j2ptr, tptr, jtstart, jchars, ichars, y, x, zchars, wp)
        if chk == false
            throw(error("data didn't check out!"))
        end

        # construct tmpvars
        dtv = DrillingTmpVars(maxj1length(j1ptr), model, Float64)

        return new{M,ETV,ITup,XT}(model, j1ptr, j2ptr, tptr, jtstart, jchars, ichars, y, x, zchars, dtv)
    end
end

DataDrill(d::AbstractDataDrill) = _data(d)
DataDrill(g::AbstractDataStructure) = DataDrill(_data(g))
reward(x::DataDrill) = reward(_model(x))

function DataDrill(d::AbstractDataDrill, m::AbstractDrillModel)
    statespace(d) == statespace(m) || throw(error("state space not same"))
    return DataDrill(
        m, j1ptr(d), j2ptr(d), tptr(d), jtstart(d),
        j1chars(d), ichars(d),_y(d), _x(d), zchars(d)
    )
end

@deprecate DataDrill(m::AbstractDrillModel, d::AbstractDataDrill) DataDrill(d,m)

# What is an observation?
#------------------------------------------

struct ObservationDrill{M<:AbstractDrillModel,ITup<:Tuple,ZTup<:Tuple,XT<:Number} <: AbstractObservation
    model::M
    ichars::ITup
    z::ZTup
    y::Int
    x::XT
end

const ObservationDynamicDrill = ObservationDrill{<:AbstractDynamicDrillModel}
const ObservationStaticDrill = ObservationDrill{ <:AbstractStaticDrillModel}

# ichars
@inline ichars(obs::ObservationDrill) = obs.ichars
# @deprecate _ichars(obs::ObservationDrill) ichars(obs)
@inline geology(obs::ObservationDrill{<:AbstractDrillModel,<:NTuple{2,Real}}) = first(ichars(obs))
@inline royalty(obs::ObservationDrill{<:AbstractDrillModel,<:NTuple{2,Real}}) = last(ichars(obs))

# zchars
zchars(obs::ObservationDrill) = obs.z
@deprecate _z(     obs::ObservationDrill) zchars(obs)
@inline logprice(  obs::ObservationDrill) = logprice(  zchars(obs))
@inline logrigrate(obs::ObservationDrill) = logrigrate(zchars(obs))
@inline year(      obs::ObservationDrill) = year(      zchars(obs))
@inline price(     obs::ObservationDrill) = exp(logprice(obs))
@inline rigrate(   obs::ObservationDrill) = exp(logrigrate(obs))

function Observation(d::AbstractDataDrill, i::Integer, j::Integer, t::Integer)
    0 < i <= length(d) || throw(BoundsError())
    j in j1_range(d,i) || j == j2ptr(d,i) || throw(BoundsError())
    t in trange(d,j)   || throw(BoundsError())
    zt = t - tstart(d,j) + jtstart(d,j)
    return ObservationDrill(_model(d), ichars(d,i), zchars(d,zt), _y(d,t), _x(d,t))
end

Observation(d::AbstractDataDrill, i, r::AbstractRegimeType, j, t) = Observation(d,i,j,t)

# @deprecate action(obs::ObservationDrill) _y(obs)
# @deprecate state(obs::ObservationDrill) _x(obs)

# API for DataDrill
#------------------------------------------

# functions in case we want to work w/out DataDrill
hasj1ptr(   j1ptr::Vector{Int}) = length(j1ptr) > 0
maxj1length(j1ptr::Vector{Int}) = hasj1ptr(j1ptr) ? maximum( diff(j1ptr) ) : 1

# access DataDrill fields
j1ptr(   d::AbstractDataDrill) = d.j1ptr
j2ptr(   d::AbstractDataDrill) = d.j2ptr
tptr(    d::AbstractDataDrill) = d.tptr
zchars(  d::AbstractDataDrill) = d.zchars
jtstart( d::AbstractDataDrill) = d.jtstart
ichars(  d::AbstractDataDrill) = d.ichars
j1chars( d::AbstractDataDrill) = d.j1chars
DrillingTmpVars(d::DataDrill) = d.dtv
DrillingTmpVars(data) = DrillingTmpVars(_data(data))

statespace(d::DataDrillPrimitive) = d.wp
statespace(d::DataDrill) = statespace(_model(d))

coefnames(x::DataDrill)  = coefnames(_model(x))

# length of things
hasj1ptr(   d::AbstractDataDrill) = hasj1ptr(j1ptr(d))
length(     d::AbstractDataDrill) = length(ichars(d))
num_i(d::AbstractDataDrill) = length(d)
maxj1length(d::AbstractDataDrill) = maxj1length(j1ptr(d))
_nparm(     d::AbstractDataDrill) = _nparm(_model(d))
j1length(   d::AbstractDataDrill) = length(j1chars(d))
@deprecate total_leases(d) j1length(d)

# getindex in fields of AbstractDataDrill
ichars( d::AbstractDataDrill, i) = getindex(ichars(d),  i)
j1ptr(  d::AbstractDataDrill, i) = getindex(j1ptr(d),   i)
j2ptr(  d::AbstractDataDrill, i) = getindex(j2ptr(d),   i)

tptr(   d::AbstractDataDrill, j) = getindex(tptr(d),    j)
jtstart(d::AbstractDataDrill, j) = getindex(jtstart(d), j)
j1chars(d::AbstractDataDrill, j) = getindex(j1chars(d), j)

_y(     d::AbstractDataDrill, t) = getindex(_y(d), t)
_x(     d::AbstractDataDrill, t) = getindex(_x(d), t)
zchars( d::AbstractDataDrill, t) = getindex(zchars(d), t)

# Iteration through j1
j1start( d::AbstractDataDrill, i) = j1ptr(d,i)
j1stop(  d::AbstractDataDrill, i) = j1start(d,i+1)-1
j1length(d::AbstractDataDrill, i) = hasj1ptr(d) ? j1stop(d,i) - j1start(d,i) + 1 : 0
j1_range(d::AbstractDataDrill, i) = hasj1ptr(d) ? (j1start(d,i) : j1stop(d,i)) : (1:0)

# iteration through t
tstart( d::AbstractDataDrill, j) = tptr(d,j)
tstop(  d::AbstractDataDrill, j) = tstart(d,j+1)-1
trange( d::AbstractDataDrill, j) = tstart(d,j) : tstop(d,j)
tlength(d::AbstractDataDrill, j) = tstop(d,j) - tstart(d,j) + 1

# iteration through zchars
zcharsvec(d::AbstractDataDrill, t0) = view(zchars(d), t0:length(zchars(d)))

# @deprecate j2_index(     data::DataDrill, i::Integer) j2ptr(  data,i)
# @deprecate j1_indexrange(data::DataDrill, i::Integer) j1range(data,i)
# @deprecate tend(         data::DataDrill, j::Integer) tstop(  data,j)
# @deprecate ilength(      data::DataDrill)             length( data)

# Unit (first layer of iteration)
#------------------------------------------

const DrillUnit = ObservationGroup{<:AbstractDataDrill}

j1length( g::DrillUnit) = j1length(_data(g), _i(g))
j1_range( g::DrillUnit) = j1_range(_data(g), _i(g))
j1start(  g::DrillUnit) = j1start( _data(g), _i(g))
j1stop(   g::DrillUnit) = j1stop(  _data(g), _i(g))
j2ptr(    g::DrillUnit) = j2ptr(   _data(g), _i(g))
ichars(   g::DrillUnit) = ichars(  _data(g), _i(g))
j1chars(  g::DrillUnit) = view(j1chars(_data(g)), j1_range(g))
uniti(    g::DrillUnit) = _i(g)

j1_sample(g::DrillUnit) = first(sample(j1_range(g), 1))

firstindex(grp::DrillUnit) = InitialDrilling()
lastindex( grp::DrillUnit) = DevelopmentDrilling()
length(    grp::DrillUnit) = DevelopmentDrilling()
eachindex( grp::DrillUnit) = (InitialDrilling(), DevelopmentDrilling())

# Convenience Constructors
InitialDrilling(    d::DrillUnit) = ObservationGroup(d,InitialDrilling())
DevelopmentDrilling(d::DrillUnit) = ObservationGroup(d,DevelopmentDrilling())

num_initial_leases(d::DrillUnit) = j1length(d)

function max_state(grp::DrillUnit)
    d = _data(grp)
    _t1range = tstart(d, j1start(grp)) : tstop(d, j1stop(grp))
    _t2range = trange(d, j2ptr(grp))
    x1 = length(_t1range) > 0 ? maximum(view(_x(d), _t1range)) : 0
    x2 = length(_t2range) > 0 ? maximum(view(_x(d), _t2range)) : 0
    return max(x1,x2)
end

max_states(d::DataDrill) = [max_state(g) for g in d]

total_wells_drilled(d::DataDrill) = [_D(statespace(_model(d)), max_state(g)) for g in d]

# Regime (second layer of iteration)
#------------------------------------------

# At the Unit level
const AbstractDrillRegime = ObservationGroup{<:DrillUnit}
const DrillInitial        = ObservationGroup{<:DrillUnit,InitialDrilling}
const DrillDevelopment    = ObservationGroup{<:DrillUnit,DevelopmentDrilling}

length(    g::DrillInitial) = j1length(_data(g))
eachindex( g::DrillInitial) = j1_range(_data(g))
firstindex(g::DrillInitial) = j1start( _data(g))
lastindex( g::DrillInitial) = j1stop(  _data(g))
j1chars(   g::DrillInitial) = j1chars( _data(g))

length(    g::DrillDevelopment) = 1
eachindex( g::DrillDevelopment) = j2ptr(_data(g))
firstindex(g::DrillDevelopment) = j2ptr(_data(g))
lastindex( g::DrillDevelopment) = j2ptr(_data(g))
j1chars(   g::DrillDevelopment) = 1

uniti(g::ObservationGroup) = uniti(_data(g))

# Lease (third layer of iteration)
#------------------------------------------

const DrillLease = ObservationGroup{<:AbstractDrillRegime}

length(    g::DrillLease) = tlength(DataDrill(g), _i(g))
eachindex( g::DrillLease) = trange( DataDrill(g), _i(g))
tstart(    g::DrillLease) = tstart( DataDrill(g), _i(g))
firstindex(g::DrillLease) = tstart(g)
lastindex( g::DrillLease) = tstop(  DataDrill(g), _i(g))

_y(g::DrillLease) = view(_y(DataDrill(g)), eachindex(g))
_x(g::DrillLease) = view(_x(DataDrill(g)), eachindex(g))

jtstart(g::DrillLease) = jtstart(DataDrill(g), _i(g))
ztrange(g::DrillLease) = (jtstart(g)-1) .+ OneTo(length(g))
zchars( g::DrillLease) = view( zchars(DataDrill(g)), ztrange(g))
ichars( g::DrillLease) = ichars(DataDrill(g), uniti(g))
j1chars(g::DrillLease) = j1chars(DataDrill(g), _i(g))
_j(g::DrillLease) = _i(g)
_regime(g::DrillLease) = _i(_data(g))

lease_ever_drilled(g::DrillLease) = sum(_y(g)) > 0

function lease_expired(g::DrillLease)
    m = _model(_data(g))
    wp = statespace(m)
    return exploratory_terminal(wp) in _x(g)
end

function getindex(g::DrillLease, t)
    Observation(DataDrill(g), uniti(g), _regime(g) ,_j(g), t)
end

# Simulation of data
# -----------------------------------------------------

is_development(lease::ObservationGroup{<:DrillInitial}) = false
is_development(lease::ObservationGroup{<:DrillDevelopment}) = true

function randsumtoone(n)
    x = rand(n)
    x ./= sum(x)
    return x
end

function initialize_x!(x, m, lease)
    x .= (2*is_development(lease)-1) .* abs.(x)
end
update_x!(x, t, m, state, d) = nothing

# check that VF is not all zeros...
check_vf_not_zero(m::AbstractStaticDrillModel) = true
check_vf_not_zero(m::AbstractDynamicDrillModel) = any(EV(value_function(m)) .!= 0)

function simulate_lease(lease::DrillLease, theta::AbstractVector{<:Number}, sim::SimulationDraw)
    m = _model(DataDrill(lease))
    length(theta) == _nparm(m) || throw(DimensionMismatch())

    tmpgrad = similar(theta)

    nper = length(lease)
    if nper > 0

        zc = zchars(lease)
        ic = ichars(lease)
        x = _x(lease)
        y = _y(lease)

        i = uniti(lease)
        ubv = Vector{Float64}(undef, length(actionspace(m)))

        initialize_x!(x, m, lease)

        dograd = false
        update_interpolation!(value_function(m), dograd)

        for t in 1:nper
            obs = ObservationDrill(m, ic, zc[t], y[t], x[t])
            f(d) = full_payoff!(tmpgrad, d, obs, theta, sim, dograd)
            actions = actionspace(obs)
            resize!(ubv, length(actions))
            ubv .= f.(actions)
            logsumexp!(ubv)
            cumsum!(ubv, ubv)
            choice = searchsortedfirst(ubv, rand())-1
            y[t] = choice
            update_x!(x, t, m, x[t], choice)
        end
    end
end

# ExogTimeVarsSample(m::AbstractDrillModel, nt::Integer) = throw(error("not defined for $(m)"))
function ExogTimeVarsSample(m::AbstractDrillModel, nt::Integer)
    timespan = range(Date(2003,10); step=Quarter(1), length=nt)
    timevars = [(x,) for x in randn(nt)]
    return ExogTimeVars(timevars, timespan)
end

# ichars_sample(m::AbstractDrillModel, num_i) = throw(error("not defined for $(m)"))
function ichars_sample(m::AbstractDrillModel, num_i)
    [(x,) for x in sample(0:1, num_i)]
end

xsample(d::UnivariateDistribution, nobs::Integer) = rand(d, nobs)
xsample(d::UnitRange, nobs::Integer) = sample(d, nobs)

function DataDrill(u::Vector, v::Vector, _zchars::ExogTimeVars, _ichars::Vector{<:NTuple{N,Number}},
    m::AbstractDrillModel, theta::AbstractVector;
    minmaxleases::UnitRange=0:3, nper_initial::UnitRange=1:10,
    nper_development::UnitRange=0:10,
    tstart::UnitRange=5:15,
    xdomain::D=Normal()
) where {D,N}

    all(u .!= v) || throw(error("u,v must be different!"))

    num_i = length(u)
    num_i == length(v) || throw(DimensionMismatch())
    num_zt = length(_zchars)

    # initial leases per unit
    initial_leases_per_unit = sample(minmaxleases, num_i)
    _jchars = vcat(collect(randsumtoone(lpu) for lpu in initial_leases_per_unit )...)
    num_initial_leases = length(_jchars)

    # observations per lease
    obs_per_lease = vcat(
        sample(nper_initial, num_initial_leases),
        sample(nper_development, num_i)
    )

    # pointers to observations
    _tptr  = 1 .+ cumsum(vcat(0, obs_per_lease))
    _j1ptr = 1 .+ cumsum(vcat(0, initial_leases_per_unit))
    _j2ptr = (last(_j1ptr)-1) .+ (1:num_i)

    nobs = last(_tptr)-1
    x = xsample(xdomain, nobs)
    y = zeros(Int, nobs)
    _jtstart = sample(tstart, num_initial_leases + num_i)

    choice_set = actionspace(m)

    data = DataDrill(m, _j1ptr, _j2ptr, _tptr, _jtstart, _jchars, _ichars, y, x, _zchars)

    θρ = theta_drill_ρ(reward(m),theta)

    # update leases
    println("Simulating $num_i units.")
    @showprogress 1 for (i,unit) in enumerate(data)
        sim = SimulationDraw(u[i], v[i], θρ)
        solve_vf_and_update_itp!(m, theta, ichars(unit), false)
        check_vf_not_zero(m) || @warn "Value Function is all 0s!"
        for regimes in unit
            for lease in regimes
                simulate_lease(lease, theta, sim)
            end
        end
    end

    return data
end

function DataDrill(m, theta; num_i=100, kwargs...)
    u = randn(num_i)
    v = randn(num_i)
    DataDrill(u,v,m,theta; kwargs...)
end

function DataDrill(u, v, m, theta; num_zt=30, kwargs...)
    _zchars = ExogTimeVarsSample(m, num_zt)
    num_i = length(u)
    _ichars = ichars_sample(m,num_i)
    return DataDrill(u,v,_zchars, _ichars,m,theta;kwargs...)
end

# ---------------------------
# DataDrlll - start only
# ---------------------------

function DataDrillStartOnly(d::AbstractDataDrill)
    _tptr = zero(tptr(d))
    _tptr[1] = 1

    for unit in d
        for regime in unit
            for lease in regime
                _tptr[_i(lease)+1] = length(lease) > 0
            end
        end
    end
    cumsum!(_tptr, _tptr)

    XT = eltype(_x(d))
    x = Vector{XT}(undef, last(_tptr)-1)
    y = Vector{Int}(undef, length(x))

    for unit in d
        for regime in unit
            for lease in regime
                if length(lease) > 0
                    trng = range_i_to_ip1(_tptr, _i(lease))
                    x[trng] .= first(_x(lease))
                    y[trng] .= first(_y(lease))
                end
            end
        end
    end

    return DataDrill(
        _model(d),
        j1ptr(d), j2ptr(d), _tptr, jtstart(d), j1chars(d),
        ichars(d), y, x, zchars(d)
    )
end

# create new DataDrill with different reward & statespace
function DataDrill(d::DataDrill, rwrd::AbstractStaticPayoff, wp::AbstractStateSpace)

    # create new DDM from old
    ddm = DynamicDrillModel(_model(d), rwrd, wp)

    if wp == statespace(d)
        newx = _x(d)
    else
        states = map(x -> state(statespace(d), x), _x(d))
        newx = map(s -> state_idx(wp, s), states)
    end

    dnew = DataDrill(
        ddm,
        j1ptr(d), j2ptr(d), tptr(d), jtstart(d), j1chars(d),
        ichars(d), _y(d), newx, zchars(d)
    )
    return dnew
end


text:
function DiffusionDefinition.diffusion_message(::Val{:LotkaVolterra})
    println("* * * * * * * * * * *")
    println("Lotka–Volterra system")
    println("* * * * * * * * * * *")
    println("Description")
    println("------------")
    println("A simple, scalar-valued predator-prey model.")
    println("It is represented by a two-dimensional process (Xₜ, Yₜ).")
    println("dXₜ = (αXₜ - βXₜYₜ)dt + σ₁dW¹ₜ,")
    println("dYₜ = (δXₜYₜ - γYₜ)dt + σ₂dW²ₜ.")
    println("The model has six parameters:")
    println("\tα, β, γ, δ, σ1, σ2")
    println("\nConstructors")
    println("------------")
    println("You may initialize the process either by passing the parameters")
    println("in the order above, as in:")
    println("α, β, γ, δ, σ1, σ2 = ...")
    println("P = LotkaVolterra(α, β, γ, δ, σ1, σ2)")
    println("or by keyword arguments")
    println("p = (α=..., β=..., γ=..., δ=..., σ1=..., σ2=...)")
    println("P = LotkaVolterra(;p...)")
    println()
end

@diffusion_process LotkaVolterra{K} begin
    :dimensions
    process --> 2
    wiener --> 2

    :parameters
    (α, β, γ, δ, σ1, σ2) --> K

    #=
    :conjugate
    phi(t, x) --> (
        (0.0, 0.0),
        (x[1], 0.0),
        (-x[1]*x[2], 0.0),
        (0.0, x[1]*x[2]),
        (0.0, -x[2]),
        (0.0, 0.0),
        (0.0, 0.0)
    )
    nonhypo(x) --> x
    num_non_hypo --> 2
    =#

    :additional
    constdiff --> true
    statespace --> LowerBoundedStateSpace((1,2),(0.0, 0.0))
end

function DiffusionDefinition.b(t, x, P::LotkaVolterra)
    @SVector [
        P.α*x[1] - P.β*x[1]*x[2],
        P.δ*x[1]*x[2] - P.γ*x[2],
    ]
end

DiffusionDefinition.σ(t, x, P::LotkaVolterra) = SDiagonal(P.σ1, P.σ2)


text:
a = 999
b = 999
for i = 999:-1:100
    s = a*b
    if string(s) == reverse(string(s))
        println(s)
        break
    else
    end
    global a = i
    if string(s) == reverse(string(s))
        println(s)
        break
    else
    end
    global b = i
    if string(s) == reverse(string(s))
        println(s)
        break
    else
    end
end


text:
using Test
import ParallelStencil
using ParallelStencil.ParallelKernel
import ParallelStencil.ParallelKernel: @reset_parallel_kernel, @is_initialized, SUPPORTED_PACKAGES, PKG_CUDA, PKG_THREADS
import ParallelStencil.ParallelKernel: @require, @prettystring, @gorgeousstring
import ParallelStencil.ParallelKernel: checkargs_parallel, checkargs_parallel_indices, parallel_indices
using ParallelStencil.ParallelKernel.Exceptions
TEST_PACKAGES = SUPPORTED_PACKAGES
@static if PKG_CUDA in TEST_PACKAGES
    import CUDA
    if !CUDA.functional() TEST_PACKAGES = filter!(x->x≠PKG_CUDA, TEST_PACKAGES) end
end

@static for package in TEST_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. parallel macros" begin
            @require !@is_initialized()
            @init_parallel_kernel($package, Float64)
            @require @is_initialized()
            @testset "@parallel" begin
                @static if $package == $PKG_CUDA
                    call = @prettystring(1, @parallel f(A))
                    @test occursin("CUDA.@cuda blocks = ParallelStencil.ParallelKernel.compute_nblocks(length.(ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A))), ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A))))) threads = ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))) f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[3])))", call)
                    @test occursin("CUDA.synchronize()", call)
                    call = @prettystring(1, @parallel ranges f(A))
                    @test occursin("CUDA.@cuda blocks = ParallelStencil.ParallelKernel.compute_nblocks(length.(ParallelStencil.ParallelKernel.promote_ranges(ranges)), ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ranges)))) threads = ParallelStencil.ParallelKernel.compute_nthreads(length.(ParallelStencil.ParallelKernel.promote_ranges(ranges))) f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[3])))", call)
                    call = @prettystring(1, @parallel nblocks nthreads f(A))
                    @test occursin("CUDA.@cuda blocks = nblocks threads = nthreads f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[3])))", call)
                    call = @prettystring(1, @parallel ranges nblocks nthreads f(A))
                    @test occursin("CUDA.@cuda blocks = nblocks threads = nthreads f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[3])))", call)
                    call = @prettystring(1, @parallel nblocks nthreads stream=mystream f(A))
                    @test occursin("CUDA.@cuda blocks = nblocks threads = nthreads stream = mystream f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[3])))", call)
                elseif $package == $PKG_THREADS
                    @test @prettystring(1, @parallel f(A)) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[3])))"
                    @test @prettystring(1, @parallel ranges f(A)) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[3])))"
                    @test @prettystring(1, @parallel nblocks nthreads f(A)) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.compute_ranges(nblocks .* nthreads)))[3])))"
                    @test @prettystring(1, @parallel ranges nblocks nthreads f(A)) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ranges), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ranges))[3])))"
                    @test @prettystring(1, @parallel stream=mystream f(A)) == "f(A, ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[1])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[2])), (Int64)(length((ParallelStencil.ParallelKernel.promote_ranges(ParallelStencil.ParallelKernel.get_ranges(A)))[3])))"
                end;
            end;
            @testset "@parallel_indices" begin
                @testset "addition of range arguments" begin
                    expansion = @gorgeousstring(1, @parallel_indices (ix,iy) f(a::T, b::T) where T <: Union{Array{Float32}, Array{Float64}} = (println("a=$a, b=$b)"); return))
                    @test occursin("f(a::T, b::T, ranges::Tuple{UnitRange, UnitRange, UnitRange}, rangelength_x::Int64, rangelength_y::Int64, rangelength_z::Int64", expansion)
                end
                @testset "Data.Array to Data.DeviceArray" begin
                    @static if $package == $PKG_CUDA
                            expansion = @prettystring(1, @parallel_indices (ix,iy) f(A::Data.Array, B::Data.Array, c::T) where T <: Integer = (A[ix,iy] = B[ix,iy]^c; return))
                            @test occursin("f(A::Data.DeviceArray, B::Data.DeviceArray,", expansion)
                    end
                end
                @testset "@parallel_indices (1D)" begin
                    A  = @zeros(4)
                    @parallel_indices (ix) function write_indices!(A)
                        A[ix] = ix;
                        return
                    end
                    @parallel write_indices!(A);
                    @test all(Array(A) .== [ix for ix=1:size(A,1)])
                end;
                @testset "@parallel_indices (2D)" begin
                    A  = @zeros(4, 5)
                    @parallel_indices (ix,iy) function write_indices!(A)
                        A[ix,iy] = ix + (iy-1)*size(A,1);
                        return
                    end
                    @parallel write_indices!(A);
                    @test all(Array(A) .== [ix + (iy-1)*size(A,1) for ix=1:size(A,1), iy=1:size(A,2)])
                end;
                @testset "@parallel_indices (3D)" begin
                    A  = @zeros(4, 5, 6)
                    @parallel_indices (ix,iy,iz) function write_indices!(A)
                        A[ix,iy,iz] = ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2);
                        return
                    end
                    @parallel write_indices!(A);
                    @test all(Array(A) .== [ix + (iy-1)*size(A,1) + (iz-1)*size(A,1)*size(A,2) for ix=1:size(A,1), iy=1:size(A,2), iz=1:size(A,3)])
                end;
                @testset "@parallel_indices (1D in 3D)" begin
                    A  = @zeros(4, 5, 6)
                    @parallel_indices (iy) function write_indices!(A)
                        A[1,iy,1] = iy;
                        return
                    end
                    @parallel 1:size(A,2) write_indices!(A);
                    @test all(Array(A)[1,:,1] .== [iy for iy=1:size(A,2)])
                end;
                @testset "@parallel_indices (2D in 3D)" begin
                    A  = @zeros(4, 5, 6)
                    @parallel_indices (ix,iz) function write_indices!(A)
                        A[ix,end,iz] = ix + (iz-1)*size(A,1);
                        return
                    end
                    @parallel (1:size(A,1), 1:size(A,3)) write_indices!(A);
                    @test all(Array(A)[:,end,:] .== [ix + (iz-1)*size(A,1) for ix=1:size(A,1), iz=1:size(A,3)])
                end;
            end;
            @testset "@parallel_async" begin
                @static if $package == $PKG_CUDA
                    call = @prettystring(1, @parallel_async f(A))
                    @test !occursin("CUDA.synchronize()", call)
                end;
            end;
            @testset "@synchronize" begin
                @static if $package == $PKG_CUDA
                    @test @prettystring(1, @synchronize()) == "CUDA.synchronize()"
                end;
            end;
            @reset_parallel_kernel()
        end;
        @testset "2. parallel macros (literal conversion)" begin
            @testset "@parallel_indices (Float64)" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, Float64)
                @require @is_initialized()
                expansion = @prettystring(@parallel_indices (ix) f!(A) = (A[ix] = A[ix] + 1.0f0; return))
                @test occursin("A[ix] = A[ix] + 1.0\n", expansion)
                @reset_parallel_kernel()
            end;
            @testset "@parallel_indices (Float32)" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, Float32)
                @require @is_initialized()
                expansion = @prettystring(@parallel_indices (ix) f!(A) = (A[ix] = A[ix] + 1.0; return))
                @test occursin("A[ix] = A[ix] + 1.0f0\n", expansion)
                @reset_parallel_kernel()
            end;
            @testset "@parallel_indices (Float16)" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, Float16)
                @require @is_initialized()
                expansion = @prettystring(@parallel_indices (ix) f!(A) = (A[ix] = A[ix] + 1.0; return))
                @test occursin("A[ix] = A[ix] + Float16(1.0)\n", expansion)
                @reset_parallel_kernel()
            end;
            @testset "@parallel_indices (ComplexF64)" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, ComplexF64)
                @require @is_initialized()
                expansion = @prettystring(@parallel_indices (ix) f!(A) = (A[ix] = 2.0f0 - 1.0f0im - A[ix] + 1.0f0; return))
                @test occursin("A[ix] = ((2.0 - 1.0im) - A[ix]) + 1.0\n", expansion)
                @reset_parallel_kernel()
            end;
            @testset "@parallel_indices (ComplexF32)" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, ComplexF32)
                @require @is_initialized()
                expansion = @prettystring(@parallel_indices (ix) f!(A) = (A[ix] = 2.0 - 1.0im - A[ix] + 1.0; return))
                @test occursin("A[ix] = ((2.0f0 - 1.0f0im) - A[ix]) + 1.0f0\n", expansion)
                @reset_parallel_kernel()
            end;
            @testset "@parallel_indices (ComplexF16)" begin
                @require !@is_initialized()
                @init_parallel_kernel($package, ComplexF16)
                @require @is_initialized()
                expansion = @prettystring(@parallel_indices (ix) f!(A) = (A[ix] = 2.0 - 1.0im - A[ix] + 1.0; return))
                @test occursin("A[ix] = ((Float16(2.0) - Float16(1.0) * im) - A[ix]) + Float16(1.0)\n", expansion)
                @reset_parallel_kernel()
            end;
        end
        @testset "3. parallel macros (numbertype ommited)" begin
            @require !@is_initialized()
            @init_parallel_kernel(package = $package)
            @require @is_initialized
            @testset "Data.Array{T} to Data.DeviceArray{T}" begin
                @static if $package == $PKG_CUDA
                        expansion = @prettystring(1, @parallel_indices (ix,iy) f(A::Data.Array{T}, B::Data.Array{T}, c<:Integer) where T <: Union{Float32, Float64}  = (A[ix,iy] = B[ix,iy]^c; return))
                        @test occursin("f(A::Data.DeviceArray{T}, B::Data.DeviceArray{T},", expansion)
                end
            end
            @reset_parallel_kernel()
        end
        @testset "4. Exceptions" begin
            @require !@is_initialized()
            @init_parallel_kernel($package, Float64)
            @require @is_initialized
            @testset "arguments @parallel" begin
                @test_throws ArgumentError checkargs_parallel();                                                    # Error: isempty(args)
                @test_throws ArgumentError checkargs_parallel(:(f()), :(something));                                # Error: last arg is not function call.
                @test_throws ArgumentError checkargs_parallel(:(f()=99));                                           # Error: last arg is not function call.
                #TODO: Here I am: kw for calls look very different: head :kw - fix in parallel.jl/shared.jl
                @test_throws ArgumentError checkargs_parallel(:(f(;s=1)));                                          # Error: function call with keyword argument.
                @test_throws ArgumentError checkargs_parallel(:ranges, :nblocks, :nthreads, :something, :(f()));    # Error: length(posargs) > 3
                @test_throws KeywordArgumentError checkargs_parallel(:(blocks=blocks), :(f()));                     # Error: blocks keyword argument is not allowed
                @test_throws KeywordArgumentError checkargs_parallel(:(threads=threads), :(f()));                   # Error: threads keyword argument is not allowed
            end;
            @testset "arguments @parallel_indices" begin
                @test_throws ArgumentError checkargs_parallel_indices();                                            # Error: length(args) != 2
                @test_throws ArgumentError checkargs_parallel_indices(:(f()=99));                                   # Error: length(args) != 2
                @test_throws ArgumentError checkargs_parallel_indices(:((ix,iy,iz)), :(f()=99), :(something));      # Error: length(args) != 2
                @test_throws ArgumentError checkargs_parallel_indices(:ix, :iy, :iz, :(f()=99));                    # Error: length(args) != 2
                @test_throws ArgumentError checkargs_parallel_indices(:(f()=99), :((ix,iy,iz)));                    # Error: last arg is not function.
                @test_throws ArgumentError checkargs_parallel_indices(:((ix,iy,iz)), :(f()));                       # Error: last arg is not function.
                @test_throws ArgumentError checkargs_parallel_indices(:((ix,iy,iz)), :(f(;s=1)=(99*s; return)))     # Error: function defines keyword.
                @test_throws ArgumentError parallel_indices(:((ix,iy,iz)), :(f()=99))                               # Error: no return statement in function.
                @test_throws ArgumentError parallel_indices(:((ix,iy,iz)), :(f()=(99; return something)))           # Error: function does not return nothing.
                #TODO: this tests does not pass anymore for unknown reasons:
                #@test_throws ArgumentError parallel_indices(:((ix,iy,iz)), :(f()=(99; if x return y end; return)))  # Error: function contains more than one return statement.
            end;
            @reset_parallel_kernel()
        end;
    end;
)) end == nothing || true;


text:
# AUTO GENERATED FILE - DO NOT EDIT

export ''_dashecharts

"""
    ''_dashecharts(;kwargs...)

A DashECharts component.

Keyword arguments:
- `id` (String; optional): The ID used to identify this component in Dash callbacks.
- `bmap_token` (String; optional)
- `brush_data` (Dict; optional)
- `click_data` (Dict; optional)
- `event` (Dict; optional)
- `fun_effects` (Array; optional)
- `fun_keys` (Array; optional)
- `fun_paths` (Dict; optional)
- `fun_prepares` (Array; optional)
- `fun_values` (Array; optional)
- `funs` (Dict; optional)
- `mapbox_token` (String; optional)
- `maps` (Dict; optional)
- `n_clicks` (Real; optional)
- `n_clicks_timestamp` (Real; optional)
- `option` (Dict; optional)
- `reset_id` (Real; optional)
- `resize_id` (Real; optional)
- `selected_data` (Dict; optional)
- `style` (Dict; optional)
"""
function ''_dashecharts(; kwargs...)
        available_props = Symbol[:id, :bmap_token, :brush_data, :click_data, :event, :fun_effects, :fun_keys, :fun_paths, :fun_prepares, :fun_values, :funs, :mapbox_token, :maps, :n_clicks, :n_clicks_timestamp, :option, :reset_id, :resize_id, :selected_data, :style]
        wild_props = Symbol[]
        return Component("''_dashecharts", "DashECharts", "dash_echarts", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Covarians/dash-echarts

max_issues_repo_name: Covarians/dash-echarts

max_forks_repo_name: Covarians/dash-echarts