text:
"""
    HDirProdPs(x,y)

Calculates horizontal direct product of two matrices with equal number of rows.
              z[i,:] is the Kronecker product of x[i,:] and y[i,:]

# Usage
z = HDirProdPs(x,y)

# Input
- `x::Array`:      T x Kx matrix
- `y::Array`:      T x Ky matrix

# Output
- `z::Array`:      T x (Kx*Ky)

# Example
```julia
julia>  x = [ 1 2;
              3 4]
julia>  y = [5 6 1;
             7 8 1]
julia> HDirProdPs(x,y)
2x6 Array{Int64,2}:
  5   6  1  10  12  2
 21  24  3  28  32  4
```

"""
function HDirProdPs(x,y)

  Kx = size(x,2)       #columns in x
  Ky = size(y,2)       #columns in y

  z = repeat(y,1,Kx) .* kron(x,ones(Int,1,Ky))  #Int: more general, small perf penalty

  return z

end
#----------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode


text:
# ------------------------------------------------------------------------------
# finite element bases
# ------------------------------------------------------------------------------

"""
Finite element basis for function spaces and test spaces
"""
abstract type AbstractBasis end

# ------------------------------------------------------------------------------
# basis types
# ------------------------------------------------------------------------------

"""
```julia
TensorBasis(
    numbernodes1d,
    numberquadraturepoints1d,
    numbercomponents,
    dimension,
    nodes1d,
    quadraturepoints1d,
    quadratureweights1d,
    interpolation1d,
    gradient1d
)
```

Tensor product basis

# Arguments:
- `numbernodes1d`:             number of nodes in 1 dimension
- `numberquadraturepoints1d`:  number of quadrature points in 1 dimension
- `numbercomponents`:          number of components
- `dimension`:                 dimension of the basis
- `nodes1d`:                   coordinates of the nodes in 1 dimension
- `quadraturepoints1d`:        coordinates of the quadrature points in 1
                                   dimension
- `quadratureweights1d`:       quadrature weights in 1 dimension
- `interpolation1d`:           interpolation matrix from nodes to quadrature
                                   points in 1 dimension
- `gradient1d`:                gradient matrix from nodes to quadrature points in
                                   1 dimension

# Returns:
- Tensor product basis object
"""
mutable struct TensorBasis <: AbstractBasis
    # data never changed
    numbernodes1d::Int
    numberquadratuepoints1d::Int
    numbercomponents::Int
    dimension::Int
    nodes1d::AbstractArray{Float64,1}
    quadraturepoints1d::AbstractArray{Float64,1}
    quadratureweights1d::AbstractArray{Float64,1}
    interpolation1d::AbstractArray{Float64,2}
    gradient1d::AbstractArray{Float64,2}
    volume::Float64
    numberelements1d::Int

    # data empty until assembled
    nodes::AbstractArray{Float64}
    quadraturepoints::AbstractArray{Float64}
    quadratureweights::AbstractArray{Float64,1}
    interpolation::AbstractArray{Float64,2}
    gradient::AbstractArray{Float64,2}
    numbermodes::Int
    modemap::AbstractArray{Int,1}
    primalnodes::AbstractArray{Int,1}
    interfacenodes::AbstractArray{Int,1}

    # inner constructor
    TensorBasis(
        numbernodes1d::Int,
        numberquadraturepoints1d::Int,
        numbercomponents::Int,
        dimension::Int,
        nodes1d::AbstractArray{Float64,1},
        quadraturepoints1d::AbstractArray{Float64,1},
        quadratureweights1d::AbstractArray{Float64,1},
        interpolation1d::AbstractArray{Float64,2},
        gradient1d::AbstractArray{Float64,2};
        numberelements1d::Int = 1,
    ) = (
        # validity checking
        if numbernodes1d < 1
            error("number of basis nodes in 1d must be at least 1") # COV_EXCL_LINE
        end;
        if numberquadraturepoints1d < 1
            error("number of quadrature points in 1d must be at least 1") # COV_EXCL_LINE
        end;
        if numbercomponents < 1
            error("number of components must be at least 1") # COV_EXCL_LINE
        end;
        if dimension < 1
            error("dimension must be at least 1") # COV_EXCL_LINE
        end;
        if length(nodes1d) != numbernodes1d
            error("must include numbernodes1d nodes") # COV_EXCL_LINE
        end;
        if length(quadraturepoints1d) != numberquadraturepoints1d
            error("must include numberquadraturepoints1d quadrature points") # COV_EXCL_LINE
        end;
        if length(quadratureweights1d) != numberquadraturepoints1d
            error("must include numberquadraturepoints1d quadrature weights") # COV_EXCL_LINE
        end;
        if size(interpolation1d) != (numberquadraturepoints1d, numbernodes1d)
            # COV_EXCL_START
            error(
                "interpolation matrix must have dimensions (numberquadraturepoints1d, numbernodes1d)",
            )
            # COV_EXCL_STOP
        end;
        if size(gradient1d) != (numberquadraturepoints1d, numbernodes1d)
            # COV_EXCL_START
            error(
                "gradient matrix must have dimensions (numberquadraturepoints1d, numbernodes1d)",
            ) # COV_EXCL_STOP
        end;

        # constructor
        new(
            numbernodes1d,
            numberquadraturepoints1d,
            numbercomponents,
            dimension,
            nodes1d,
            quadraturepoints1d,
            quadratureweights1d,
            interpolation1d,
            gradient1d,
            (maximum(nodes1d) - minimum(nodes1d))^dimension,
            numberelements1d,
        )
    )
end

# printing
# COV_EXCL_START
function Base.show(io::IO, basis::TensorBasis)
    print(
        io,
        basis.numberelements1d == 1 ? "" : "macro-element ",
        "tensor product basis:\n    numbernodes1d: ",
        basis.numbernodes1d,
        "\n    numberquadraturepoints1d: ",
        basis.numberquadratuepoints1d,
        "\n    numbercomponents: ",
        basis.numbercomponents,
    )
    if basis.numberelements1d != 1
        print(io, "\n    numberelements1d: ", basis.numberelements1d)
    end
    print(io, "\n    dimension: ", basis.dimension)
end
# COV_EXCL_STOP

"""
```julia
NonTensorBasis(
    numbernodes,
    numberquadraturepoints,
    numbercomponents,
    dimension,
    nodes,
    quadraturepoints,
    quadratureweights,
    interpolation,
    gradient
)
```
Non-tensor basis

# Arguments:
- `numbernodes`:             number of nodes
- `numberquadraturepoints`:  number of quadrature points
- `numbercomponents`:        number of components
- `dimension`:               dimension of the basis
- `nodes`:                   coordinates of the nodes
- `quadraturepoints`:        coordinates of the quadrature points
- `quadratureweights`:       quadrature weights
- `interpolation`:           interpolation matrix from nodes to quadrature points
- `gradient`:                gradient matrix from nodes to quadrature points

# Returns:
- Non-tensor product basis object
"""
mutable struct NonTensorBasis <: AbstractBasis
    # data never changed
    numbernodes::Int
    numberquadraturepoints::Int
    numbercomponents::Int
    dimension::Int
    nodes::AbstractArray{Float64}
    quadraturepoints::AbstractArray{Float64,1}
    quadratureweights::AbstractArray{Float64,1}
    interpolation::AbstractArray{Float64,2}
    gradient::AbstractArray{Float64,2}
    volume::Float64
    numbermodes::Int
    modemap::AbstractArray{Int,1}
    numberelements::Int
    primalnodes::AbstractArray{Int,1}
    interfacenodes::AbstractArray{Int,1}

    # inner constructor
    NonTensorBasis(
        numbernodes::Int,
        numberquadraturepoints::Int,
        numbercomponents::Int,
        dimension::Int,
        nodes::Int,
        quadraturepoints::AbstractArray{Float64},
        quadratureweights::AbstractArray{Float64,1},
        interpolation::AbstractArray{Float64,2},
        gradient::AbstractArray{Float64,2};
        numberelements::Int = 1,
    ) = (
        # validity checking
        if numbernodes < 1
            error("number of nodes must be at least 1") # COV_EXCL_LINE
        end;
        if numberquadraturepoints < 1
            error("number of quadrature points must be at least 1") # COV_EXCL_LINE
        end;
        if numbercomponents < 1
            error("number of components must be at least 1") # COV_EXCL_LINE
        end;
        if dimension < 1
            error("dimension must be at least 1") # COV_EXCL_LINE
        end;
        if size(nodes) != (numbernodes, dimension)
            error("must include (numbernodes, dimension) nodal coordinates") # COV_EXCL_LINE
        end;
        if length(quadraturepoints) != (numberquadraturepoints, dimension)
            error("must include (numberquadraturepoints, dimension) quadrature points") # COV_EXCL_LINE
        end;
        if length(quadratureweights) != numberquadraturepoints
            error("must include sufficient quadrature weights") # COV_EXCL_LINE
        end;
        if size(interpolation) != (numberquadraturepoints, numbernodes)
            # COV_EXCL_START
            error(
                "interpolation matrix must have dimensions (numberquadraturepoints, numbernodes)",
            )
            # COV_EXCL_STOP
        end;
        if size(gradient) != (q * dimension, numbernodes)
            # COV_EXCL_START
            error(
                "gradient matrix must have dimensions (numberquadraturepoints*dimension, numbernodes)",
            )
            # COV_EXCL_STOP
        end;
        if length(modemap) != numbercomponents * numbernodes
            error("must map the modes for each basis node") # COV_EXCL_LINE
        end;

        # compute volume
        volume = 1;
        for d = 1:dimension
            volume *= (maximum(nodes[:, d]) - minimum(nodes[:, d]))
        end;

        # constructor
        new(
            numbernodes,
            numberquadraturepoints,
            numbercomponents,
            dimension,
            nodes,
            quadraturepoints,
            quadratureweights,
            interpolation,
            gradient,
            volume,
            maximum(modemap),
            modemap,
            numberelements,
        )
    )
end

# printing
# COV_EXCL_START
function Base.show(io::IO, basis::NonTensorBasis)
    print(
        io,
        basis.numberelements == 1 ? "" : "macro-element ",
        "non-tensor product basis:\n    numbernodes: ",
        basis.numbernodes,
        "\n    numberquadraturepoints: ",
        basis.numberquadratuepoints,
        "\n    numbercomponents: ",
        basis.numbercomponents,
    )
    if basis.numberelements != 1
        print(io, "\n    numberelements: ", basis.numberelements)
    end
    print(io, "\n    dimension: ", basis.dimension)
end
# COV_EXCL_STOP

# ------------------------------------------------------------------------------
# basis properties
# ------------------------------------------------------------------------------

"""
```julia
getnumbernodes(basis)
```

Get the number of nodes for the basis

# Arguments:
- `basis`:  basis to compute number of nodes

# Returns:
- Integer number of basis nodes

# Example:
```jldoctest
# get number of nodes for basis
basis = TensorH1LagrangeBasis(4, 3, 2, 2);

# note: either syntax works
numbernodes = LFAToolkit.getnumbernodes(basis);
numbernodes = basis.numbernodes;

# verify
@assert numbernodes == 4^2

# output

```
"""
function getnumbernodes(basis::TensorBasis)
    return basis.numbernodes1d^basis.dimension
end

"""
```julia
getnodes(basis)
```

Get nodes for basis

# Returns:
- Basis nodes array

# Arguments:
- `basis`:  basis to compute nodes

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get basis quadrature weights
    basis = TensorH1LagrangeBasis(4, 3, 2, dimension);

    # note: either syntax works
    nodes = LFAToolkit.getnodes(basis);
    nodes = basis.nodes;

    # verify
    truenodes1d = [-1, -√(1/5), √(1/5), 1];
    truenodes = [];
    if dimension == 1
        truenodes = truenodes1d;
    elseif dimension == 2
        truenodes =
            transpose(hcat([[[x, y] for x in truenodes1d, y in truenodes1d]...]...));
    elseif dimension == 3
        truenodes = transpose(hcat([[
            [x, y, z] for x in truenodes1d, y in truenodes1d, z in truenodes1d
        ]...]...));
    end

    @assert truenodes ≈ nodes
end

# output

```
"""
function getnodes(basis::TensorBasis)
    # assembled if needed
    if !isdefined(basis, :nodes)
        nodes = []
        if basis.dimension == 1
            # 1D
            nodes = basis.nodes1d
        elseif basis.dimension == 2
            # 2D
            nodes =
                transpose(hcat([[[x, y] for x in basis.nodes1d, y in basis.nodes1d]...]...))
        elseif basis.dimension == 3
            # 3D
            nodes = transpose(
                hcat(
                    [
                        [
                            [x, y, z] for x in basis.nodes1d, y in basis.nodes1d,
                            z in basis.nodes1d
                        ]...,
                    ]...,
                ),
            )
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        basis.nodes = nodes
    end

    # return
    return getfield(basis, :nodes)
end

"""
```julia
getnumberquadraturepoints(basis)
```

Get the number of quadrature points for the basis

# Arguments:
- `basis`:  basis to compute number of quadrature points

# Returns:
- Integer number of basis quadrature points

# Example:
```jldoctest
# get number of quadrature points for basis
basis = TensorH1LagrangeBasis(4, 3, 2, 2);

# note: either syntax works
numberquadraturepoints = LFAToolkit.getnumberquadraturepoints(basis);
numberquadraturepoints = basis.numberquadraturepoints;

# verify
@assert numberquadraturepoints == 3^2

# output

```
"""
function getnumberquadraturepoints(basis::TensorBasis)
    return basis.numberquadratuepoints1d^basis.dimension
end

"""
```julia
getquadraturepoints(basis)
```

Get quadrature points for basis

# Returns:
- Basis quadrature points array

# Arguments:
- `basis`: basis to compute quadrature points

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get basis quadrature weights
    basis = TensorH1LagrangeBasis(4, 3, 2, dimension);

    # note: either syntax works
    quadraturepoints = LFAToolkit.getquadraturepoints(basis);
    quadraturepoints = basis.quadraturepoints;

    # verify
    truequadraturepoints1d = [-√(3/5), 0, √(3/5)];
    truequadraturepoints = [];
    if dimension == 1
        truequadraturepoints = truequadraturepoints1d;
    elseif dimension == 2
        truequadraturepoints = transpose(hcat([[
            [x, y] for x in truequadraturepoints1d, y in truequadraturepoints1d
        ]...]...));
    elseif dimension == 3
        truequadraturepoints = transpose(hcat([[
            [x, y, z]
            for
            x in truequadraturepoints1d,
            y in truequadraturepoints1d, z in truequadraturepoints1d
        ]...]...));
    end

    @assert truequadraturepoints ≈ quadraturepoints
end

# output

```
"""
function getquadraturepoints(basis::TensorBasis)
    # assembled if needed
    if !isdefined(basis, :quadraturepoints)
        quadraturepoints = []
        if basis.dimension == 1
            # 1D
            quadraturepoints = basis.quadraturepoints1d
        elseif basis.dimension == 2
            # 2D
            quadraturepoints = transpose(
                hcat(
                    [
                        [
                            [x, y] for x in basis.quadraturepoints1d,
                            y in basis.quadraturepoints1d
                        ]...,
                    ]...,
                ),
            )
        elseif basis.dimension == 3
            # 3D
            quadraturepoints = transpose(
                hcat(
                    [
                        [
                            [x, y, z] for x in basis.quadraturepoints1d,
                            y in basis.quadraturepoints1d, z in basis.quadraturepoints1d
                        ]...,
                    ]...,
                ),
            )
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        basis.quadraturepoints = quadraturepoints
    end

    # return
    return getfield(basis, :quadraturepoints)
end

"""
```julia
getquadratureweights(basis)
```

Get full quadrature weights vector for basis

# Returns:
- Basis quadrature weights vector

# Arguments:
- `basis`:  basis to compute quadrature weights

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get basis quadrature weights
    basis = TensorH1LagrangeBasis(4, 3, 2, dimension);

    # note: either syntax works
    quadratureweights = LFAToolkit.getquadratureweights(basis);
    quadratureweights = basis.quadratureweights;

    # verify
    trueweights1d = [5/9, 8/9, 5/9];
    trueweights = [];
    if dimension == 1
        trueweights = trueweights1d;
    elseif dimension == 2
        trueweights = kron(trueweights1d, trueweights1d);
    elseif dimension == 3
        trueweights = kron(trueweights1d, trueweights1d, trueweights1d);
    end

    @assert trueweights ≈ quadratureweights
end

# output

```
"""
function getquadratureweights(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :quadratureweights)
        quadratureweights = []
        if basis.dimension == 1
            # 1D
            quadratureweights = basis.quadratureweights1d
        elseif basis.dimension == 2
            # 2D
            quadratureweights = kron(basis.quadratureweights1d, basis.quadratureweights1d)
        elseif basis.dimension == 3
            # 3D
            quadratureweights = kron(
                basis.quadratureweights1d,
                basis.quadratureweights1d,
                basis.quadratureweights1d,
            )
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        basis.quadratureweights = quadratureweights
    end

    # return
    return getfield(basis, :quadratureweights)
end

"""
```julia
getinterpolation(basis)
```

Get full interpolation matrix for basis

# Arguments:
- `basis`:  basis to compute interpolation matrix

# Returns:
- Basis interpolation matrix

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get basis interpolation matrix
    basis = TensorH1LagrangeBasis(4, 3, 2, dimension);

    # note: either syntax works
    interpolation = LFAToolkit.getinterpolation(basis);
    interpolation = basis.interpolation;

    # verify
    for i in 1:3^dimension
        total = sum(interpolation[i, :]);
        @assert total ≈ 1.0
    end
end

# output

```
"""
function getinterpolation(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :interpolation)
        interpolation = []
        if basis.dimension == 1
            # 1D
            interpolation = kron(I(basis.numbercomponents), basis.interpolation1d)
        elseif basis.dimension == 2
            # 2D
            interpolation = kron(
                I(basis.numbercomponents),
                basis.interpolation1d,
                basis.interpolation1d,
            )
        elseif basis.dimension == 3
            # 3D
            interpolation = kron(
                I(basis.numbercomponents),
                basis.interpolation1d,
                basis.interpolation1d,
                basis.interpolation1d,
            )
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        basis.interpolation = interpolation
    end

    # return
    return getfield(basis, :interpolation)
end

"""
```julia
getgradient(basis)
```

Get full gradient matrix for basis

# Arguments:
- `basis`:  basis to compute gradient matrix

# Returns:
- Basis gradient matrix

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get basis gradient matrix
    basis = TensorH1LagrangeBasis(4, 3, 2, dimension);

    # note: either syntax works
    gradient = LFAToolkit.getgradient(basis);
    gradient = basis.gradient;

    # verify
    for i in 1:dimension*3^dimension
        total = sum(gradient[i, :]);
        @assert abs(total) < 1e-14
    end
end

# output

```
"""
function getgradient(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :gradient)
        gradient = []
        if basis.dimension == 1
            # 1D
            gradient = kron(I(basis.numbercomponents), basis.gradient1d)
        elseif basis.dimension == 2
            # 2D
            gradient = kron(
                I(basis.numbercomponents),
                [
                    kron(basis.gradient1d, basis.interpolation1d)
                    kron(basis.interpolation1d, basis.gradient1d)
                ],
            )
        elseif basis.dimension == 3
            # 3D
            gradient = kron(
                I(basis.numbercomponents),
                [
                    kron(basis.gradient1d, basis.interpolation1d, basis.interpolation1d)
                    kron(basis.interpolation1d, basis.gradient1d, basis.interpolation1d)
                    kron(basis.interpolation1d, basis.interpolation1d, basis.gradient1d)
                ],
            )
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        basis.gradient = gradient
    end

    # return
    return getfield(basis, :gradient)
end

"""
```julia
getnumbermodes(basis)
```

Get number of modes for basis

# Arguments:
- `basis`:  basis to compute number of modes

# Returns:
- Number of modes for basis

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get number of basis modes
    basis = TensorH1LagrangeBasis(4, 3, 2, dimension);

    # note: either syntax works
    numbermodes = LFAToolkit.getnumbermodes(basis);
    numbermodes = basis.numbermodes;

    # verify
    @assert numbermodes == 2*3^dimension
end

# output

```
"""
function getnumbermodes(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :modemap)
        basis.numbermodes = maximum(basis.modemap)
    end

    # return
    return getfield(basis, :numbermodes)
end

"""
```julia
getmodemap(basis)
```

Get mode mapping vector for basis

# Arguments:
- `basis`:  basis to compute mode map vector

# Returns:
- Basis mode map vector

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get mode map vector
    basis = TensorH1LagrangeBasis(4, 3, 1, dimension);

    # note: either syntax works
    modemap = LFAToolkit.getmodemap(basis);
    modemap = basis.modemap;

    # verify
    truemodemap1d = [1, 2, 3, 1];
    truemodemap = [];
    if dimension == 1
        truemodemap = truemodemap1d;
    elseif dimension == 2
        truemodemap = [[
            i + (j - 1)*3 for i in truemodemap1d, j in truemodemap1d
        ]...];
    elseif dimension == 3
        truemodemap = [[
            i +
            (j - 1)*3 +
            (k - 1)*3^2
            for i in truemodemap1d, j in truemodemap1d, k in truemodemap1d
        ]...];
    end

    @assert truemodemap == modemap
end

# output

```
"""
function getmodemap(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :modemap)
        modemap1d = [1:basis.numbernodes1d;]
        modemap1d[end] = 1
        modemap = []
        if basis.dimension == 1
            # 1D
            modemap = modemap1d
        elseif basis.dimension == 2
            # 2D
            modemap = [
                [
                    i + (j - 1) * (basis.numbernodes1d - 1) for i in modemap1d,
                    j in modemap1d
                ]...,
            ]
        elseif basis.dimension == 3
            # 3D
            modemap = [
                [
                    i +
                    (j - 1) * (basis.numbernodes1d - 1) +
                    (k - 1) * (basis.numbernodes1d - 1)^2 for i in modemap1d,
                    j in modemap1d, k in modemap1d
                ]...,
            ]
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        numbermodes1component = maximum(modemap)
        modemap = vcat(
            [
                ((i - 1) * numbermodes1component) .+ modemap for
                i = 1:basis.numbercomponents
            ]...,
        )
        basis.modemap = modemap
        basis.numbermodes = maximum(modemap)
    end

    # return
    return getfield(basis, :modemap)
end

"""
```julia
getprimalnodes(basis)
```

Get primal nodes for basis

# Arguments:
- `basis`:  basis to compute primal nodes

# Returns:
- Basis primal nodes vector

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get mode map vector
    p = 4
    basis = TensorH1LagrangeBasis(p, p, 1, dimension);

    # note: either syntax works
    primalnodes = LFAToolkit.getprimalnodes(basis);
    primalnodes = basis.primalnodes;

    # verify
    trueprimalnodes = []
    if dimension == 1
        trueprimalnodes = [1, p];
    elseif dimension == 2
        trueprimalnodes = [1, p, p^2-p+1, p^2]
    elseif dimension == 3
        trueprimalnodes = [1, p, p^2-p+1, p^2, p^3-p^2+1, p^3-p^2+p, p^3-p+1, p^3]
    end

    @assert trueprimalnodes == primalnodes
end

# output

```
"""
function getprimalnodes(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :primalnodes)
        primalnodes1component = []
        p = basis.numbernodes1d
        if basis.dimension == 1
            # 1D
            primalnodes1component = [1, p]
        elseif basis.dimension == 2
            # 2D
            primalnodes1component = [1, p, p^2 - p + 1, p^2]
        elseif basis.dimension == 3
            # 3D
            primalnodes1component =
                [1, p, p^2 - p + 1, p^2, p^3 - p^2 + 1, p^3 - p^2 + p, p^3 - p + 1, p^3]
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        primalnodes = vcat(
            [
                (i - 1) * basis.numbernodes .+ primalnodes1component for
                i = 1:basis.numbercomponents
            ]...,
        )
        basis.primalnodes = primalnodes
    end

    # return
    return getfield(basis, :primalnodes)
end

"""
```julia
getinterfacenodes(basis)
```

Get interface nodes for basis

# Arguments:
- `basis`:  basis to compute primal nodes

# Returns:
- Basis primal nodes vector

# Example:
```jldoctest
# test for all supported dimensions
for dimension in 1:3
    # get mode map vector
    p = 4
    basis = TensorH1LagrangeBasis(p, p, 1, dimension);

    # note: either syntax works
    interfacenodes = LFAToolkit.getinterfacenodes(basis);
    interfacenodes = basis.interfacenodes;

    # verify
    trueinterfacenodes = []
    if dimension == 1
        trueinterfacenodes = [1, p];
    elseif dimension == 2
        trueinterfacenodes = [1:p..., p^2-p+1:p^2...]
        for i = 1:p-2
            push!(trueinterfacenodes, i*p+1)
            push!(trueinterfacenodes, (i+1)*p)
        end
    elseif dimension == 3
        trueinterfacenodes = [1:p^2..., p^3-p^2+1:p^3...]
        for i = 1:p-2
            push!(trueinterfacenodes, i*p^2+1:i*p^2+p...)
            push!(trueinterfacenodes, (i+1)*p^2-p+1:(i+1)*p^2...)
            push!(trueinterfacenodes, i*p^2+p+1:p:(i+1)*p^2-2*p+1...)
            push!(trueinterfacenodes, i*p^2+2*p:p:(i+1)*p^2-p...)
        end
    end
    trueinterfacenodes = sort(trueinterfacenodes)

    @assert trueinterfacenodes == interfacenodes
end

# output

```
"""
function getinterfacenodes(basis::TensorBasis)
    # assemble if needed
    if !isdefined(basis, :interfacenodes)
        interfacenodes1component = []
        p = basis.numbernodes1d
        if basis.dimension == 1
            # 1D
            interfacenodes1component = [1, p]
        elseif basis.dimension == 2
            # 2D
            interfacenodes1component = [1:p..., p^2-p+1:p^2...]
            for i = 1:p-2
                push!(interfacenodes1component, i * p + 1)
                push!(interfacenodes1component, (i + 1) * p)
            end
        elseif basis.dimension == 3
            # 3D
            interfacenodes1component = [1:p^2..., p^3-p^2+1:p^3...]
            for i = 1:p-2
                push!(interfacenodes1component, i*p^2+1:i*p^2+p...)
                push!(interfacenodes1component, (i+1)*p^2-p+1:(i+1)*p^2...)
                push!(interfacenodes1component, i*p^2+p+1:p:(i+1)*p^2-2*p+1...)
                push!(interfacenodes1component, i*p^2+2*p:p:(i+1)*p^2-p...)
            end
        else
            throw(DomainError(basis.dimension, "Dimension must be less than or equal to 3")) # COV_EXCL_LINE
        end
        interfacenodes1component = sort(interfacenodes1component)
        interfacenodes = vcat(
            [
                (i - 1) * basis.numbernodes .+ interfacenodes1component for
                i = 1:basis.numbercomponents
            ]...,
        )
        basis.interfacenodes = interfacenodes
    end

    # return
    return getfield(basis, :interfacenodes)
end

"""
```julia
getnumberelements(basis)
```

Get the number of elements for the basis

# Arguments:
- `basis`:  basis to compute number of micro-elements

# Returns:
- Integer number of basis micro-elements

# Example:
```jldoctest
# get number of nodes for basis
basis = TensorH1LagrangeMacroBasis(4, 4, 1, 2, 2);

# note: either syntax works
numbernodes = LFAToolkit.getnumberelements(basis);
numbernodes = basis.numberelements;

# verify
@assert numbernodes == 2^2

# output

```
"""
function getnumberelements(basis::TensorBasis)
    return basis.numberelements1d^basis.dimension
end

# ------------------------------------------------------------------------------
# get/set property
# ------------------------------------------------------------------------------

function Base.getproperty(basis::TensorBasis, f::Symbol)
    if f == :numbernodes
        return getnumbernodes(basis)
    elseif f == :numberquadraturepoints
        return getnumberquadraturepoints(basis)
    elseif f == :nodes
        return getnodes(basis)
    elseif f == :quadratureweights
        return getquadratureweights(basis)
    elseif f == :quadraturepoints
        return getquadraturepoints(basis)
    elseif f == :interpolation
        return getinterpolation(basis)
    elseif f == :gradient
        return getgradient(basis)
    elseif f == :numbermodes
        return getnumbermodes(basis)
    elseif f == :modemap
        return getmodemap(basis)
    elseif f == :primalnodes
        return getprimalnodes(basis)
    elseif f == :interfacenodes
        return getinterfacenodes(basis)
    elseif f == :numberelements
        return getnumberelements(basis)
    else
        return getfield(basis, f)
    end
end

function Base.setproperty!(basis::TensorBasis, f::Symbol, value)
    if f == :numbernodes1d ||
       f == :numberquadraturepoints1d ||
       f == :numbercomponents ||
       f == :dimension ||
       f == :nodes1d ||
       f == :quadraturepoints1d ||
       f == :quadratureweights1d ||
       f == :interpolation1d ||
       f == :gradient1d ||
       f == :volume
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(basis, f, value)
    end
end

function Base.setproperty!(basis::NonTensorBasis, f::Symbol, value)
    if f == :numbernodes ||
       f == :numberquadraturepoints ||
       f == :numbercomponents ||
       f == :dimension ||
       f == :nodes ||
       f == :quadraturepoints ||
       f == :quadratureweights ||
       f == :interpolation ||
       f == :gradient ||
       f == :volume ||
       f == :numbermodes ||
       f == :modemap
        throw(ReadOnlyMemoryError()) # COV_EXCL_LINE
    else
        return setfield!(basis, f, value)
    end
end

# ------------------------------------------------------------------------------
# gradient on stretched mesh
# ------------------------------------------------------------------------------

"""
```julia
getdXdxgradient(basis, mesh)
```

Get gradient adjusted for mesh stretching

# Arguments:
- `basis`:  basis to compute gradient
- `mesh`:   mesh to compute gradient

# Returns:
- gradient matrix multiplied by change of coordinates adjoint

# Example:
```jldoctest
for dimension = 1:3
    # mesh
    mesh = []
    if dimension == 1
        mesh = Mesh1D(2.0)
    elseif dimension == 2
        mesh = Mesh2D(2.0, 3.0)
    elseif dimension == 3
        mesh = Mesh3D(2.0, 3.0, 4.0)
    end

    # basis
    basis = TensorH1LagrangeBasis(4, 3, 1, dimension);

    # get gradient on mesh
    gradient = LFAToolkit.getdXdxgradient(basis, mesh);

    # verify
    nodes = basis.nodes;
    linearfunction = [];
    truegradient = [];
    if dimension == 1
        linearfunction = nodes;
        truegradient = [1/2*ones(basis.numberquadraturepoints)...]
    elseif dimension == 2
        linearfunction = (nodes[:, 1] + nodes[:, 2]);
        truegradient = [
            1/2*ones(basis.numberquadraturepoints)...
            1/3*ones(basis.numberquadraturepoints)...
        ]
    elseif dimension == 3
        linearfunction = (nodes[:, 1] + nodes[:, 2] + nodes[:, 3]);
        truegradient = [
            1/2*ones(basis.numberquadraturepoints)...
            1/3*ones(basis.numberquadraturepoints)...
            1/4*ones(basis.numberquadraturepoints)...
        ]
    end

    @assert gradient*linearfunction ≈ truegradient
end

# output

```
"""
function getdXdxgradient(basis::TensorBasis, mesh::Mesh)
    # check compatibility
    dimension = basis.dimension
    if (dimension == 1 && typeof(mesh) != Mesh1D) ||
       (dimension == 2 && typeof(mesh) != Mesh2D) ||
       (dimension == 3 && typeof(mesh) != Mesh3D)
        error("mesh dimension must match basis dimension") # COV_EXCL_LINE
    end

    # coordinate transformation
    gradient = basis.gradient

    # adjust for mesh
    if dimension == 1
        # 1D
        return gradient / mesh.dx
    elseif dimension == 2
        # 2D
        numberquadraturepoints = basis.numberquadraturepoints
        return [
            gradient[1:numberquadraturepoints, :] / mesh.dx
            gradient[numberquadraturepoints+1:end, :] / mesh.dy
        ]
    elseif dimension == 3
        # 3D
        numberquadraturepoints = basis.numberquadraturepoints
        return [
            gradient[1:numberquadraturepoints, :] / mesh.dx
            gradient[numberquadraturepoints+1:2*numberquadraturepoints, :] / mesh.dy
            gradient[2*numberquadraturepoints+1:end, :] / mesh.dz
        ]
    end
end

function getdXdxgradient(basis::NonTensorBasis, mesh::Mesh)
    throw(error("dXdxgradient unimplemented for non-tensor bases"))
end

"""
```julia
getdxdXquadratureweights(basis, mesh)
```

Get quadrature weights adjusted for mesh stretching

# Arguments:
- `basis`:  basis to compute quadratureweights
- `mesh`:   mesh to compute quadratureweights

# Returns:
- quadrature weights multiplied by change of coordinates adjoint

# Example:
```jldoctest
mesh = Mesh1D(2.0)

# basis
basis = TensorH1LagrangeBasis(4, 3, 1, 1);

# get gradient on mesh
weights = LFAToolkit.getdxdXquadratureweights(basis, mesh);

# verify
@assert basis.quadratureweights * mesh.volume / basis.volume ≈ weights

# output

```
"""
function getdxdXquadratureweights(basis::TensorBasis, mesh::Mesh)
    # check compatibility
    dimension = basis.dimension
    if (dimension == 1 && typeof(mesh) != Mesh1D) ||
       (dimension == 2 && typeof(mesh) != Mesh2D) ||
       (dimension == 3 && typeof(mesh) != Mesh3D)
        error("mesh dimension must match basis dimension") # COV_EXCL_LINE
    end

    # adjust for mesh
    return basis.quadratureweights * mesh.volume / basis.volume
end

function getdxdXquadratureweights(basis::NonTensorBasis, mesh::Mesh)
    throw(error("dxdXquadratureweights unimplemented for non-tensor bases"))
end

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

ext: jl


text:
using PyPlot, Revise, transD_GP, Random, Statistics
## from matlab circle code
matlabLM = vec([5.588e-9  3.9533e-9  2.6558e-9  1.7276e-9  1.1292e-9  7.3569e-10  4.9102e-10  3.401e-10  2.4425e-10  1.8129e-10  1.3582e-10  1.0205e-10  7.7079e-11  5.8448e-11  4.4735e-11  3.4357e-11 2.6101e-11  1.9288e-11])
matlabHM = vec([5.4669e-10  3.4682e-10  2.4341e-10  1.7922e-10  1.3526e-10  1.0387e-10  8.0645e-11  6.343e-11  5.0158e-11  3.9315e-11  3.0087e-11  2.2234e-11  1.5781e-11  1.0723e-11  6.9601e-12 4.3111e-12  2.5472e-12  1.4351e-12  7.7144e-13  3.9605e-13  2.0261e-13])
## model
zfixed   = [-1e5,   0,    20,   50]
rho      = [1e12,   10,   1,   100]
nmax = 100
##  geometry and modeling parameters
rRx = 17.
zRx = -37.0
zTx = -35.0
rTx = 12.607
lowpassfcs = [300000, 450000.0]
ntimesperdecade = 10
nfreqsperdecade = 5
# Note that the receiver depth needs to be in same model layer as transmitter.
## LM times and ramp
LM_times = [
1.539e-5  1.939e-5  2.439e-5  3.139e-5  3.939e-5  4.939e-5  6.239e-5  7.839e-5  9.939e-5  0.00012539  0.00015739  0.00019939  0.00025039  0.00031539  0.00039739  0.00050039  0.00063039  0.00079339
 1.9e-5    2.4e-5    3.1e-5    3.9e-5    4.9e-5    6.2e-5    7.8e-5    9.9e-5    0.000125  0.000157    0.000199    0.00025     0.000315    0.000397    0.0005      0.00063     0.000793    0.000999]'
LM_times = exp.(mean(log.(LM_times),dims=2))[:]
LM_ramp =[
-0.001  -0.0009146  -0.0007879  -0.0005964  0.0  4.629e-7  8.751e-7  1.354e-6  2.54e-6  3.972e-6  5.404e-6  5.721e-6  6.113e-6  6.663e-6   8.068e-6  0.00125
 0.0     0.6264      0.9132      0.9905     1.0  0.9891    0.9426    0.8545    0.6053   0.303     0.04077   0.01632   0.004419  0.0006323  0.0       0.0]'
## HM times and ramp
HM_times = [
7.539e-5  9.639e-5  0.00012239  0.00015439  0.00019639  0.00024739  0.00031239  0.00039439  0.00049739  0.00062739  0.00079039  0.00099639  0.00125539  0.00158139  0.00199139  0.00250839  0.00315839  0.00397739  0.00500839  0.00630639  0.00793939
9.6e-5    0.000122  0.000154    0.000196    0.000247    0.000312    0.000394    0.000497    0.000627    0.00079     0.000996    0.001255    0.001581    0.001991    0.002508    0.003158    0.003977    0.005008    0.006306    0.007939    0.009739]'
HM_times = exp.(mean(log.(HM_times),dims=2))[:]
HM_ramp = [
-0.01  -0.008386  -0.00638  -0.003783  0.0  3.96e-7  7.782e-7  1.212e-6  3.44e-6  1.981e-5  3.619e-5  3.664e-5  3.719e-5   3.798e-5  3.997e-5  0.01
  0.0    0.4568     0.7526    0.9204    1.0  0.9984   0.9914    0.9799    0.9175   0.4587    0.007675  0.003072  0.0008319  0.000119  0.0       0.0]'
## LM operator
Flm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      lowpassfcs = lowpassfcs,
                      times  = LM_times,
                      ramp   = LM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx)
## HM operator
Fhm = transD_GP.AEM_VMD_HMD.HFieldDHT(
                      ntimesperdecade = ntimesperdecade,
                      nfreqsperdecade = nfreqsperdecade,
                      lowpassfcs = lowpassfcs,
                      times  = HM_times,
                      ramp   = HM_ramp,
                      nmax   = nmax,
                      zTx    = zTx,
                      rRx    = rRx,
                      rTx    = rTx,
                      zRx    = zRx)
## get the fields in place
transD_GP.AEM_VMD_HMD.getfieldTD!(Flm, zfixed, rho)
transD_GP.AEM_VMD_HMD.getfieldTD!(Fhm, zfixed, rho)
## plot
μ = transD_GP.AEM_VMD_HMD.μ
figure()
s1 = subplot(211)
loglog(Flm.times,μ*abs.(Flm.dBzdt), label="lm_julia")
loglog(Fhm.times,μ*abs.(Fhm.dBzdt), label="hm_julia")
loglog(Flm.times, matlabLM, ".", label="lm_matlab")
loglog(Fhm.times, matlabHM, ".", label="hm_matlab")
ylabel("dBzdt V/Am^4")
grid(which="both")
legend()
s2 = subplot(212, sharex=s1)
semilogx(Flm.times,(μ*Flm.dBzdt-matlabLM)./matlabLM*100, label="lm")
semilogx(Fhm.times,(μ*Fhm.dBzdt-matlabHM)./matlabHM*100, label="hm")
ylabel("% difference")
xlabel("time s")
grid(which="both")
legend()


text:
# This file contains the structure and common functions for SpatialWeights
mutable struct SpatialWeights
    n::Int64
    neighs::Vector{Vector{Int64}}
    weights::Vector{Vector{Float64}}
    nneighs::Vector{Int64}
    transform::Symbol
end

function SpatialWeights(W::AbstractMatrix; standardize = true)

    n, nj = size(W)

    n == nj || throw(ArgumentError("matrix is not square"))

    neighs = copy.(fill(Int[], n))
    weights = copy.(fill(Float64[], n))
    nneighs = zeros(Int, n)

    for i in 1:n
        neighs[i] = (1:n)[W[i,:] .!= 0]
        nneighs[i] = length(neighs[i])
        if standardize
            weights[i] = ones(nneighs[i]) ./ nneighs[i]
        else
            weights[i] = W[i,neighs[i]]
        end
    end

    SpatialWeights(n, neighs, weights, nneighs, standardize ? :row : :original)

end

"""
    cardinalities(W::SpatialWeights)
Return a vector with the number of neighbors of each observation.
"""
cardinalities(W::SpatialWeights)::Vector{Int64} = W.nneighs;

"""
    neighbors(W::SpatialWeights, i::Int)
Return a vector with the neighbors of ``i``.
"""
neighbors(W::SpatialWeights, i::Int)::Vector{Int64} = W.neighs[i];

"""
    weights(W::SpatialWeights, i::Int)
Return a vector with the weights of the neighbors of ``i``.
"""
weights(W::SpatialWeights, i::Int)::Vector{Float64} = W.weights[i];

"""
    wtransform!(W::SpatialWeights, style::Symbol)
In-place transformation of the weights using the specified ``style``.

# Weights transformation
`style` can be one of the following:
- `:binary`: 1 if neighbor, 0 if not.
- `:row`: row-standardized. Each row sum equals one.
"""
function wtransform!(W::SpatialWeights, style::Symbol)

    n = W.n

    (style == :binary || style == :row) ||  throw(ArgumentError("unkown transformation $(style)"))

    if style == :binary
        # Binary transformation
        for i in 1:n
            W.weights[i] = ones(length(W.neighs[i]))
        end

    else style == :row
        # Row transform
        for i in 1:n
            W.weights[i] = ones(length(W.neighs[i])) ./ length(W.neighs[i])
        end

    end

    W.transform = style
    return nothing

end

"""
    wtransform(W::SpatialWeights, style::Symbol)
Returns a transformed copy of the weights matrix.
"""
function wtransform(W::SpatialWeights, style::Symbol)::SpatialWeights
    Wnew = deepcopy(W)
    wtransform!(Wnew, style)
    return Wnew
end

"""
    wtransformation (W::SpatialWeights)
Return the current transformation of the spatial weights.
"""
wtransformation(W::SpatialWeights)::Symbol = W.transform ;

"""
    nislands (W::SpatialWeights)
Return the number of islands in the spatial weights object.
"""
nislands(W::SpatialWeights)::Int64 = sum(W.nneighs .== 0) ;

"""
    islands (W::SpatialWeights)
Return a vector with the islands in the spatial weights object.
"""
function islands(W::SpatialWeights)::Vector{Int64}
    islandsvec = Int[]
    for i = 1:W.n
        if W.nneighs[i] == 0
            push!(islandsvec, i)
        end
    end

    return islandsvec
end

# Descriptive statistics
nobs(W::SpatialWeights)::Int64 = W.n;

mean(W::SpatialWeights)::Float64 = mean(W.nneighs);

Base.minimum(W::SpatialWeights)::Int64 = minimum(W.nneighs);

Base.maximum(W::SpatialWeights)::Int64 = maximum(W.nneighs);

median(W::SpatialWeights)::Float64 = median(W.nneighs);

# Transformation functions
function Base.Matrix(W::SpatialWeights)::Matrix{Float64}
    n = W.n

    mat = zeros(n, n)
    for i in 1:n
        mat[i, W.neighs[i]] = W.weights[i]
    end

    return mat

end

function sparse(W::SpatialWeights)::SparseMatrixCSC{Float64, Int64}
    n = W.n
    tn = sum(W.nneighs)

    # Build I, J and V vectors to build the spatial matrix
    I = zeros(Int, tn)
    J = zeros(Int, tn)
    V = zeros(tn)

    nlast = 0
    for i in 1:n
        sub = (nlast + 1):(nlast + W.nneighs[i])
        I[sub] = fill(i, W.nneighs[i])
        J[sub] = W.neighs[i]
        V[sub] = W.weights[i]

        nlast = nlast + W.nneighs[i]
    end

    return sparse(I, J, V)

end

function Base.show(io::IO, W::SpatialWeights)

    print(io, "Spatial Weights \n")
    print(io, "Observations: ", nobs(W), " \n")
    print(io, "Transformation: ", wtransformation(W), "\n")
    print(io, "Minimum nunmber of neighbors: ", minimum(W), "\n")
    print(io, "Maximum nunmber of neighbors: ", maximum(W), "\n")
    print(io, "Average number of neighbors: ", round(mean(W), digits = 4), "\n")
    print(io, "Median number of neighbors: ", median(W), "\n")
    print(io, "Islands (isloated): ", nislands(W), "\n")
    print(io, "Density: ", round(100 * sum(W.nneighs) / (W.n  * W.n), digits = 4), "% \n")

end


text:
"""
    MenuAdventures.Onto

A sub-module that allows things to be on top of other things.

```jldoctest
julia> using MenuAdventures

julia> using MenuAdventures.Testing

julia> using MenuAdventures.Onto

julia> import MenuAdventures: ever_possible

julia> @universe struct Universe <: AbstractUniverse
        end;

julia> @noun struct Room <: AbstractRoom
            already_lit::Bool = true
        end;

julia> @noun struct Person <: Noun
        end;

julia> @noun struct Thing <: Noun
        end;

julia> @noun struct Table <: Noun
        end;

julia> ever_possible(::Take, ::Reachable, ::Thing) = true;

julia> ever_possible(::PutOnto, ::Reachable, ::Table) = true;

julia> ever_possible(::GoOnto, ::Immediate, ::Table) = true;

julia> cd(joinpath(pkgdir(MenuAdventures), "test", "Onto")) do
            check_choices() do interface
                you = Person(
                    "Brandon",
                    grammatical_person = second_person,
                    indefinite_article = "",
                )
                room = Room("room")
                universe = Universe(you, interface = interface)
                universe[room, you] = Containing()
                universe[room, Table("table")] = Containing()
                universe[room, Thing("thing")] = Containing()
                universe
            end
        end
true
```
"""
module Onto

using MenuAdventures: Action, get_object, get_mover, Immediate, Inventory, Reachable, Relationship, Verb
import MenuAdventures: argument_domains, ever_possible, print_sentence, string_relationship_to, verb_for

"""
    Supporting()

A is `Supporting` B means B is on top of A.

A [`Relationship`](@ref)
"""
struct Supporting <: Relationship end

export Supporting

verb_for(::Supporting) = Verb("support")

function string_relationship_to(thing_answer, ::Supporting, parent_thing)
    buffer = IOBuffer()
    print(buffer, thing_answer.text)
    print(buffer, ' ')
    print(buffer, "on ")
    show(IOContext(buffer, :is_subject => false, :subject => get_object(thing_answer)), parent_thing)
    String(take!(buffer))
end

"""
    PutOnto()

Put something from your [`Inventory`](@ref) onto something [`Reachable`](@ref).

An [`Action`](@ref).
By default, `ever_possible(::PutOnto, ::Inventory, _) = true`, that is, you can always put something
from your inventory onto a surface.
"""
struct PutOnto <: Action end

export PutOnto

ever_possible(::PutOnto, ::Inventory, _) = true

function argument_domains(::PutOnto)
    Inventory(), Reachable()
end

function (::PutOnto)(universe, thing, parent_thing)
    universe[parent_thing, thing] = Supporting()
    return false
end

function print_sentence(io, ::PutOnto, thing_answer, parent_thing_answer)
    print(io, "Put ")
    print(io, thing_answer.text)
    print(io, " onto ")
    print(io, parent_thing_answer.text)
end

"""
    GoOnto()

Go onto something [`Immediate`](@ref).

An [`Action`](@ref).
"""
struct GoOnto <: Action end

export GoOnto

function argument_domains(::GoOnto)
    (Immediate(),)
end

function (::GoOnto)(universe, place)
    PutOnto()(universe, get_mover(universe), place)
    return false
end

function print_sentence(io, ::GoOnto, place_answer)
    print(io, "Go onto ")
    print(io, place_answer.text)
end

end


text:
## PARTE 2
cd(@__DIR__)
using DelimitedFiles
input = readdlm("input.txt");

##
apply_move(t, move::AbstractString, value) = apply_move(t, Val(Symbol(move)), value)
apply_move((x, y, aim), ::Val{:forward}, value) =
    (x + value, y + aim * value, aim)
apply_move((x, y, aim), ::Val{:up}, value) =
    (x, y, aim - value)
apply_move((x, y, aim), ::Val{:down}, value) =
    (x, y, aim + value)
apply_move(t, ::Val, value) = t

##
x, y, _ = foldl(eachrow(input); init=(0,0,0)) do cur, (move, value)
    apply_move(cur, move, value)
end
@show x * y


text:
# Simulate a Poisson point process on a rectangle
# Author: H. Paul Keeler, 2019.
# Website: hpaulkeeler.com
# Repository: github.com/hpaulkeeler/posts
# For more details, see the post:
# hpaulkeeler.com/testing-the-julia-language-with-point-process-simulations/

#Note: Need the .+ for adding a scalar to an array
#Also need . for sqrt, exp, cos, sin etc and assigning scalars to arrays
#Best to use vectors instead of 1-D matrices eg x=rand(n),  NOT x=rand(n,1).

using Distributions #for random simulations
using Plots #for plotting

#Simulation window parameters
xMin=0;xMax=1;
yMin=0;yMax=1;
xDelta=xMax-xMin;yDelta=yMax-yMin; #rectangle dimensions
areaTotal=xDelta*yDelta;

#Point process parameters
lambda=100; #intensity (ie mean density) of the Poisson process

#Simulate Poisson point process
numbPoints=rand(Poisson(areaTotal*lambda)); #Poisson number of points
xx=xDelta*rand(numbPoints).+xMin;#x coordinates of Poisson points
yy=yDelta*(rand(numbPoints)).+yMin;#y coordinates of Poisson points

#Plotting
plot1=scatter(xx, yy, xlabel ="x",ylabel ="y",leg = false);
display(plot1);

max_stars_repo_name: hpkeeler/posts

max_issues_repo_name: hpkeeler/posts

max_forks_repo_name: hpkeeler/posts


text:
module DeepGetfield

export getter, findfield, denest, @deep, @deepf

using MacroTools

"""
    follow(data, operations)
Executes the `operation`s on `data` to return the field the operations are pointing at.
"""
function follow(data,operations::Vector{Expr})
    isempty(operations) && (return data)
    f = eval(operations[1])
    return follow(Base.invokelatest(f,data), operations[2:end])
end

function follow(data,operations::Vector{<:Function})
    isempty(operations) && (return data)
    f = operations[1]
    return follow(f(data), operations[2:end])
end

"""
Evaluates the expressions and returns a Function.
"""
function ops2fun(ops)
    if @capture(ops[1], x->g_(x, s_))
        ex = :(x-> $g(x, $s))
    elseif @capture(ops[1], x->g_.(x, s_))
        ex = :(x-> $(g).(x, $s))
    end
    for i = 2:length(ops)
        if @capture(ops[i], x->g_(x, s_))
            ex = :( x->$g($ex(x), $s) )
        elseif @capture(ops[i], x->g_.(x, s_))
            ex = :( x->$g.($ex(x), $s) )
        end
    end
    f = eval(ex)
end

"""
    operations::Vector{Expr} = findfield(data, field, maxdepth=8)
Recursively (depth first) searches through the composite Type `data` in search For `field`
Returns a vector of Function expressions that If followed, leads to the field.
The result should be postprocessed by `tobroadcast`.
"""
function findfield(data, field, maxdepth=8)
    operations = Expr[]
    findfield(data,field, maxdepth, operations)
end

function findfield(data, field, maxdepth, operations::Vector{Expr})
    if length(operations) > maxdepth
        return nothing
    end
    dt = typeof(follow(data, operations))
    if dt <: AbstractArray || dt <: Tuple
        push!(operations, :(x-> getindex(x,1)))
        found = findfield(data, field, maxdepth, operations)
        found != nothing && (return operations)
    else
        names = fieldnames(dt)
        for name in names
            push!(operations, :(x-> getfield(x,$(QuoteNode(name)))))
            if name == field
                return operations
            end
            found = findfield(data,field, maxdepth, operations)
            found != nothing && (return operations)
            pop!(operations)
        end
    end
    return nothing
end

"""
    ops::Vector{Expr} = tobroadcast(operations::Vector{Expr})
Searches For `getindex` operations, removes them and broadcasts subsequent operations.
"""
function tobroadcast(operations)
    operations = copy(operations)
    broadcasts = 0
    i = 1
    while i <= length(operations)
        op = operations[i]
        match = @capture(op, x->f_(x,arg_))
        if match && f == :getindex
            broadcasts += 1
            deleteat!(operations, i)
            @capture(operations[i], x->f_(x, arg_))
            if f == :getindex # TODO: only handles two consequtive getindex
                broadcasts += 1
                deleteat!(operations, i)
                @capture(operations[i], x->f_(x, arg_))
                field = QuoteNode(arg.value)
                ex = :(y -> (x -> $(f).(x,$(field))).(y))
            else
                field = QuoteNode(arg.value)
                if broadcasts == 1
                    ex = :(x -> $(f).(x,$(field)))
                else
                    ex = :(x -> $(f).(x,$(field)))
                    for j = 2:broadcasts
                        ex = :(y -> ($(ex)).(y))
                    end
                end
            end
            operations[i] = ex
        elseif broadcasts >= 1 # getfield
            field = QuoteNode(arg.value)
            if broadcasts == 1
                ex = :(x -> $(f).(x,$(field)))
            else
                ex = :(x -> $(f).(x,$(field)))
                for j = 2:broadcasts
                    ex = :(y -> ($(ex)).(y))
                end
            end
            operations[i] = ex
        end
        i += 1
    end
    operations
end

"""
    denest(x::Array{Array...})

Flattens the structure of x from an array of arrays to a tensor. Can handle deep nestings.
```julia
x = [1, 2]
y = [x, x, x]
z = [y, y, y, y]
denest(x) == x
denest(y) == [x x x]
denest(z) == cat(3, denest(y), denest(y), denest(y), denest(y))
```
"""
function denest(x)
    dims = calcdims(x)
    x = denest_rec(x)
    x = reshape(x, reverse(dims)...)
    x
end

function denest_rec(x)
    if x isa Array{<: Array}
        x = vcat(x...)
        return denest(x)
    end
    x
end

calcdims(x) = calcdims(x[1],[length(x)])
calcdims(x, dims) = dims
calcdims(x::AbstractArray, dims) = calcdims(x[1], push!(dims, length(x)))

"""
    @deep data.field
Deep getfield

Searches through the structure `data` and returns a `field` where `field` is somewhere deep inside `data`.
Arrays are broadcasted over, but the search only looks at the first element.
"""
macro deep(ex)
    @capture(ex, data_.field_) || error("Expected an expression on the form data.field")
    quote
        Base.invokelatest(getter($(esc(data)), $(QuoteNode(field))), $(esc(data)))
    end
end

"""
    @deepf data.field
Deep getfield

Searches through the structure `data` and returns a function `data -> field` where `field` is somewhere deep inside `data`.
Arrays are broadcasted over, but the search only looks at the first element.
"""
macro deepf(ex)
    @capture(ex, data_.field_) || error("Expected an expression on the form data.field")
    quote
        getter($(esc(data)), $(QuoteNode(field)))
    end
end

struct DeepGetter{F}
    f::F
    path::Vector{Symbol}
    operandtype::DataType
end # module

Base.show(io::IO, z::DeepGetter) = print(io, "Getter: ", path...)
function Base.show(io::IO, ::MIME"text/plain", g::DeepGetter)
    print(io, "DeepGetter($(g.path[end]))\nPath: ")
    for s in g.path
        print(io, s, " ")
    end
    print(io, "\nOperates on Type $(g.operandtype)\n")
end

(g::DeepGetter)(x) = g.f(x)

"""
    getter(data, field, maxdepth=8)
Function version of `@deepf`. Allows specification of maximum search depth. See also `denest`, `@deepf`.
"""
function getter(data, args...)
    # TODO: preevaluate expressions
    ops = findfield(data, args...)
    path = fieldsym.(ops)
    path = [s isa Symbol ? s : :broadcast for s ∈ path]
    ops_b = tobroadcast(ops)
    f = ops2fun(ops_b)
    DeepGetter(f, path, typeof(data))
end

function fieldsym(op)
    last = op.args[2].args[2].args[3]
    return last isa QuoteNode ? last.value : last
end

end


text:
"""
Basic continuous multivariate parameter state type

# Constructors

## BasicContMuvParameterState{N<:Real}(value::Vector{N}, <optional arguments>)

Construct a basic continuous multivariate parameter state with some ``value``.

###  Optional arguments:

* ``monitor::Vector{Bool}=fill(false, 9)``: 9-element Boolean vector indicating which of the target-related fields are
stored by the state.
* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))``: the diagnostic values of the state.

## BasicContMuvParameterState{N<:Real}(value::Vector{N}, monitor::Vector{Symbol}, <optional arguments>)

Construct a basic continuous multivariate parameter state with some ``value`` and tracked target-related fields specified by
``monitor``.

###  Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))``: the diagnostic values of the state.

## BasicContMuvParameterState{N<:Real}(size::Integer, <optional arguments>)

Construct a basic continuous multivariate parameter state with a ``value`` of specified ``size``.

###  Optional arguments:

* ``monitor::Vector{Bool}=fill(false, 9)``: 9-element Boolean vector indicating which of the target-related fields are
stored by the state.
* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``::Type{N}=Float64``: the element type of the state value.
* ``diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))``: the diagnostic values of the state.

## BasicContMuvParameterState{N<:Real}(size::Integer, monitor::Vector{Symbol}, <optional arguments>)

Construct a basic continuous multivariate parameter state with a ``value`` of specified ``size`` and tracked target-related
fields specified by ``monitor``.

###  Optional arguments:

* ``diagnostickeys::Vector{Symbol}=Symbol[]``: the diagnostic keys of the state.
* ``::Type{N}=Float64``: the element type of the state value.
* ``diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))``: the diagnostic values of the state.

# Examples

```julia
julia> state = BasicContMuvParameterState(zeros(Float64, 2), [:logtarget])
Klara.BasicContMuvParameterState{Float64}([0.0,0.0],NaN,NaN,NaN,Float64[],Float64[],Float64[],0x0 Array{Float64,2},0x0 Array{Float64,2},0x0 Array{Float64,2},0x0x0 Array{Float64,3},0x0x0 Array{Float64,3},0x0x0 Array{Float64,3},Any[],2,Symbol[])

julia> state.value
2-element Array{Float64,1}:
 0.0
 0.0
```
"""
mutable struct BasicContMuvParameterState{N<:Real} <: ParameterState{Continuous, Multivariate}
  "Vector value of basic continuous multivariate parameter state"
  value::Vector{N}
  "Value of log-likelihood at the state's value"
  loglikelihood::N
  "Value of log-prior at the state's value"
  logprior::N
  "Value of log-target at the state's value"
  logtarget::N
  "Value of gradient of log-likelihood at the state's value"
  gradloglikelihood::Vector{N}
  "Value of gradient of log-prior at the state's value"
  gradlogprior::Vector{N}
  "Value of gradient of log-target at the state's value"
  gradlogtarget::Vector{N}
  "Value of metric tensor of log-likelihood at the state's value"
  tensorloglikelihood::Matrix{N}
  "Value of metric tensor of log-prior at the state's value"
  tensorlogprior::Matrix{N}
  "Value of metric tensor of log-target at the state's value"
  tensorlogtarget::Matrix{N}
  "Value of derivatives of metric tensor of log-likelihood at the state's value"
  dtensorloglikelihood::Array{N, 3}
  "Value of derivatives of metric tensor of log-prior at the state's value"
  dtensorlogprior::Array{N, 3}
  "Value of derivatives of metric tensor of log-target at the state's value"
  dtensorlogtarget::Array{N, 3}
  "Diagnostic values associated with the sampling of the state"
  diagnosticvalues::Vector
  "Integer representing the length of vector value of basic continuous multivariate parameter state"
  size::Integer
  "Diagnostic keys associated with the sampling of the state"
  diagnostickeys::Vector{Symbol}
  diffmethods::Union{DiffMethods, Void}
  diffstate::Union{DiffState, Void}
end

function BasicContMuvParameterState(
  value::Vector{N},
  monitor::Vector{Bool}=fill(false, 9),
  diagnostickeys::Vector{Symbol}=Symbol[],
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  diagnosticvalues::Vector=Array{Any}(length(diagnostickeys)),
) where N<:Real
  v = convert(N, NaN)

  s = length(value)

  l = Array{Integer}(9)
  for i in 1:9
    l[i] = (monitor[i] == false ? zero(Integer) : s)
  end

  diffstate = DiffState()

  if diffopts != nothing
    if diffopts.order == 1
      for (i, field) in ((1, :resultll), (2, :resultlp), (3, :resultlt))
        if diffopts.targets[i]
          setfield!(diffstate, field, DiffResults.GradientResult(value))
        end
      end
    else
      for (i, field) in ((1, :resultll), (2, :resultlp), (3, :resultlt))
        if diffopts.targets[i]
          setfield!(diffstate, field, DiffResults.HessianResult(value))
        end
      end
    end

    if diffopts.mode == :reverse
      for (i, field) in ((4, :cfggll), (5, :cfgglp), (6, :cfgglt))
        if diffopts.targets[i-3]
          setfield!(diffstate, field, ReverseDiff.GradientConfig(value))
        end
      end

      if diffopts.order == 2
        for (i, field) in ((7, :cfgtll), (8, :cfgtlp), (9, :cfgtlt))
          if diffopts.targets[i-6]
            setfield!(diffstate, field, ReverseDiff.HessianConfig(value))
          end
        end
      end
    else
      if diffopts.chunksize == 0
        for (i, field, diffmethod) in ((4, :cfggll, :closurell), (5, :cfgglp, :closurelp), (6, :cfgglt, :closurelt))
          if diffopts.targets[i-3]
            setfield!(diffstate, field, ForwardDiff.GradientConfig(getfield(diffmethods, diffmethod), value))
          end
        end

        if diffopts.order == 2
          for (i, diffmethod, diffresult, diffconfig) in (
            (7, :closurell, :resultll, :cfgtll),
            (8, :closurelp, :resultlp, :cfgtlp),
            (9, :closurelt, :resultlt, :cfgtlt)
          )
            if diffopts.targets[i-6]
              setfield!(
                diffstate,
                diffconfig,
                ForwardDiff.HessianConfig(getfield(diffmethods, diffmethod), getfield(diffstate, diffresult), value)
              )
            end
          end
        end
      else
        for (i, field, diffmethod) in ((4, :cfggll, :closurell), (5, :cfgglp, :closurelp), (6, :cfgglt, :closurelt))
          if diffopts.targets[i-3]
            setfield!(
              diffstate,
              field,
              ForwardDiff.GradientConfig(
                getfield(diffmethods, diffmethod), value, ForwardDiff.Chunk{diffopts.chunksize}()
              )
            )
          end
        end

        if diffopts.order == 2
          for (i, diffmethod, diffresult, diffconfig) in (
            (7, :closurell, :resultll, :cfgtll),
            (8, :closurelp, :resultlp, :cfgtlp),
            (9, :closurelt, :resultlt, :cfgtlt)
          )
            if diffopts.targets[i-6]
              setfield!(
                diffstate,
                diffconfig,
                ForwardDiff.HessianConfig(
                  getfield(diffmethods, diffmethod),
                  getfield(diffstate, diffresult),
                  value,
                  ForwardDiff.Chunk{diffopts.chunksize}()
                )
              )
            end
          end
        end
      end
    end
  end

  BasicContMuvParameterState{N}(
    value,
    v,
    v,
    v,
    Array{N}(l[1]),
    Array{N}(l[2]),
    Array{N}(l[3]),
    Array{N}(l[4], l[4]),
    Array{N}(l[5], l[5]),
    Array{N}(l[6], l[6]),
    Array{N}(l[7], l[7], l[7]),
    Array{N}(l[8], l[8], l[8]),
    Array{N}(l[9], l[9], l[9]),
    diagnosticvalues,
    s,
    diagnostickeys,
    diffmethods,
    diffstate
  )
end

function BasicContMuvParameterState(
  value::Vector{N},
  monitor::Vector{Symbol},
  diagnostickeys::Vector{Symbol}=Symbol[],
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))
) where N<:Real
  fnames = fieldnames(BasicContMuvParameterState)
  BasicContMuvParameterState(
    value, [fnames[i] in monitor ? true : false for i in 5:13], diagnostickeys, diffmethods, diffopts, diagnosticvalues
  )
end

BasicContMuvParameterState(
  size::Integer,
  monitor::Vector{Bool}=fill(false, 9),
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))
) where {N<:Real} =
  BasicContMuvParameterState(Array{N}(size), monitor, diagnostickeys, diffmethods, diffopts, diagnosticvalues)

BasicContMuvParameterState(
  size::Integer,
  monitor::Vector{Symbol},
  diagnostickeys::Vector{Symbol}=Symbol[],
  ::Type{N}=Float64,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  diagnosticvalues::Vector=Array{Any}(length(diagnostickeys))
) where {N<:Real} =
  BasicContMuvParameterState(Array{N}(size), monitor, diagnostickeys, diffmethods, diffopts, diagnosticvalues)

value_support{N<:Real}(::Type{BasicContMuvParameterState{N}}) = Continuous
value_support(::BasicContMuvParameterState) = Continuous

variate_form{N<:Real}(::Type{BasicContMuvParameterState{N}}) = Multivariate
variate_form(::BasicContMuvParameterState) = Multivariate

eltype{N<:Real}(::Type{BasicContMuvParameterState{N}}) = N
eltype{N<:Real}(::BasicContMuvParameterState{N}) = N

generate_empty(
  state::BasicContMuvParameterState,
  diffmethods::Union{DiffMethods, Void}=nothing,
  diffopts::Union{DiffOptions, Void}=nothing,
  monitor::Vector{Bool}=[isempty(getfield(state, fieldnames(BasicContMuvParameterState)[i])) ? false : true for i in 5:13]
) =
  BasicContMuvParameterState(state.size, monitor, state.diagnostickeys, eltype(state), diffmethods, diffopts)


text:
function solve()
    a, b, c, d = [parse(Int, x) for x in split(readline())]
    max(a*c, a*d, b*c, b*d)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
export ConcreteDiscretePost

"""
    ConcreteDiscretePost <: DiscretePost

Textbook implementation of a discrete post operator, using concrete polyhedra
intersections.

### Fields

- `options` -- an `Options` structure that holds the algorithm-specific options

### Notes

This operator requires that the `Polyhedra` library is loaded,
because some concrete operations between polytopes are used.

Currently, we assume that source invariants, target invariants and guards are
polytopes in constraint representation resp. half-spaces.

### Algorithm

The algorithm is based on [Flowpipe-Guard Intersection for Reachability
Computations with Support Functions](http://spaceex.imag.fr/sites/default/files/frehser_adhs2012.pdf).
"""
struct ConcreteDiscretePost <: DiscretePost
    options::Options

    function ConcreteDiscretePost(𝑂::Options)
        𝑂copy = copy(𝑂)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:check_invariant_intersection], false)
        check_aliases_and_add_default_value!(𝑂.dict, 𝑂copy.dict, [:overapproximation], Hyperrectangle)
        return new(𝑂copy)
    end
end

# convenience constructor from pairs of symbols
ConcreteDiscretePost(𝑂::Pair{Symbol,<:Any}...) = ConcreteDiscretePost(Options(Dict{Symbol,Any}(𝑂)))

# default options for the LazyDiscretePost discrete post operator
ConcreteDiscretePost() = ConcreteDiscretePost(Options())

init(𝒫::ConcreteDiscretePost, 𝒮::AbstractSystem, 𝑂::Options) = init!(𝒫, 𝒮, copy(𝑂))

function init!(𝒫::ConcreteDiscretePost, 𝒮::AbstractSystem, 𝑂::Options)
    @assert isdefined(Main, :Polyhedra) "this algorithm needs the package " *
            "'Polyhedra' to be loaded"

    𝑂.dict[:n] = statedim(𝒮, 1)

    # solver-specific options (adds default values for unspecified options)
    𝑂out = validate_solver_options_and_add_default_values!(𝑂)

    return 𝑂out
end

function tube⋂inv(𝒫::ConcreteDiscretePost,
                  reach_tube::Vector{<:AbstractReachSet{<:LazySet{N}}},
                  invariant,
                  start_interval
                 ) where {N}
    # take intersection with source invariant

    Rsets = Vector{AbstractReachSet{<:LazySet{N}}}()
    for reach_set in reach_tube
        rs = set(reach_set)
        if rs isa CartesianProductArray && length(array(rs)) == 1
            # TODO workaround for lazy X0
            rs = overapproximate(rs, BoxDirections)
        end
        if 𝒫.options[:check_invariant_intersection] && isdisjoint(rs, invariant)
            break
        end
        push!(Rsets,
              substitute(reach_set, set=intersection(rs, invariant),
                         time_start=time_start(reach_set) + start_interval[1],
                         time_end=time_end(reach_set) + start_interval[2]))
    end

    return Rsets
end

function post(𝒫::ConcreteDiscretePost,
              HS::HybridSystem,
              waiting_list::Vector{Tuple{Int, <:AbstractReachSet{LazySet{N}}, Int}},
              passed_list,
              source_loc_id,
              tube⋂inv,
              jumps,
              options
             ) where {N}
    jumps += 1
    for trans in out_transitions(HS, source_loc_id)
        info("Considering transition: $trans")
        target_loc_id = target(HS, trans)
        target_loc = HS.modes[target(HS, trans)]
        target_invariant = target_loc.X
        constrained_map = resetmap(HS, trans)
        guard = stateset(constrained_map)

        # perform jumps
        post_jump = Vector{ReachSet{LazySet{N}}}()
        sizehint!(post_jump, length(tube⋂inv))
        for reach_set in tube⋂inv
            # check intersection with guard
            R⋂G = intersection(guard, set(reach_set))
            if isempty(R⋂G)
                continue
            end

            # apply assignment
            A⌜R⋂G⌟ = apply_assignment(𝒫, constrained_map, R⋂G)

            # intersect with target invariant
            A⌜R⋂G⌟⋂I = intersection(target_invariant, A⌜R⋂G⌟)
            if isempty(A⌜R⋂G⌟⋂I)
                continue
            end

            # store result
            push!(post_jump, ReachSet{LazySet{N}}(A⌜R⋂G⌟⋂I,
                                                     time_start(reach_set),
                                                     time_end(reach_set)))
        end

        postprocess(𝒫, HS, post_jump, options, waiting_list, passed_list,
            target_loc_id, jumps)
    end
end

# --- handling assignments ---

function apply_assignment(𝒫::ConcreteDiscretePost,
                          constrained_map::ConstrainedLinearMap,
                          R⋂G::LazySet;
                          kwargs...)
    return linear_map(constrained_map.A, R⋂G)
end


text:
# Load & initialize CUDA driver

# const libcuda = dlopen("libcuda")

macro cucall(fv, argtypes, args...)
    f = eval(fv)
    quote
        _curet = ccall( ($(Meta.quot(f)), CUDA_LIB), Cint, $argtypes, $(args...) )
        if _curet != 0
            throw(CuDriverError(Int(_curet)))
        end
    end
end

function initialize()
    @cucall(cuInit, (Cint,), 0)
    println("CUDA Driver Initialized")
end

initialize()

# Get driver version

function driver_version()
    a = Cint[0]
    @cucall(cuDriverGetVersion, (Ptr{Cint},), a)
    return Int(a[1])
end

const DriverVersion = driver_version()

if DriverVersion < 4000
    error("CUDA of version 4.0 or above is required.")
end

# box a variable into array

cubox{T}(x::T) = T[x]


text:
using DrWatson
@quickactivate "ChaosThroughBilliards"
include(srcdir("style.jl"))

# re = billiard_rectangle(x, y)
# si = billiard_sinai(y/4, x, y)
# st = billiard_stadium(x-y, y)

bd = Billiard(Antidot([0.0, 0.0], 1.5, false))

InteractiveDynamics.obfill(::Antidot) = RGBAf0(0,0,0,0)
InteractiveDynamics.obls(::Antidot) = nothing
InteractiveDynamics.obcolor(::Antidot) = to_color(:white)

N = 8
ps = [Particle(0, 0, 2π*(i/N)) for i in 1:N]

# interactive_billiard(bd, ps;
# backgroundcolor = :black, particle_size = 2.0, add_controls = false)

billiard_video(
    videodir("circle.mp4"), bd, deepcopy(ps);
    frames = 80, backgroundcolor = :black, colors = COLORS,
    tailwidth = 3.5, particle_size = 2, res = (1200, 1200),
    plot_particles = true,
    dt = 0.0001, speed = 150, tail = 10000, # this makes ultra fine temporal resolution
)

max_stars_repo_name: JuliaDynamics/ChaosThroughBilliards

max_issues_repo_name: JuliaDynamics/ChaosThroughBilliards

max_forks_repo_name: JuliaDynamics/ChaosThroughBilliards


text:
module Par

using Comm, GIDs, Foldable, Functor, FunRefs, Funs, Monad

import Base: eltype, remote
import Foldable: freduce
import Functor: fmap
import Monad: tycon, valtype, munit, mjoin, mbind

export par   # @par
export rpar, rcall   # @rpar, @rcall
export make_local, get_local
export unwrap
export remote   # @remote

# resc(expr) = esc(Base.localize_vars(expr, false))
# # Return an iterable consisting of a Callable and arguments,
# # corresponding to the content of args of a call Expr
# function make_callargs(expr)
#     if false && Funs.iscall(expr)
#         # Optimize calls by not wrapping them in a lambda expression
#         Funs.getcallargs(expr)
#     else
#         {:(()->$expr)}
#     end
# end

function par(f::Callable; R::Type=eltype(f))
    r = FunRef{R}()
    @schedule r[] = f()
    r::FunRef{R}
end

# macro par(args...)
#     if length(args) == 1
#         (expr,) = args
#         resc(:(par(()->$expr)))
#     elseif length(args) == 2
#         (R, expr) = args
#         resc(:(par(R=$R, ()->$expr)))
#     else
#         error("Expected: @par [<type>] <expr>)")
#     end
# end

immutable Item
    item::Any
end

function rpar(f::Callable, ref::FunRef; R::Type=eltype(f))
    r = FunRef{R}()
    rgid = makegid(r)
    @schedule rexec(proc(ref)) do
        rpar2(rgid, f)
    end
    r::FunRef{R}
end
rpar(f::Callable, p::Integer; R::Type=eltype(f)) = rpar(R=R, f, Int(p))
function rpar(f::Callable, p::Int; R::Type=eltype(f))
    r = FunRef{R}()
    rgid = makegid(r)
    rexec(p) do
        rpar2(rgid, f)
    end
    r::FunRef{R}
end
function rpar2(rgid::GID, f::Callable)
    res = f()
    rexec(rgid.proc) do
        rpar3(rgid, Item(res))
    end
end
function rpar3(rgid::GID, res::Item)
    r = getobj(rgid)
    free(rgid)
    r[] = res.item
end

# TODO: Introduce LocalFunRef (aka Future) that does not obtain a GID
# for the item
# TODO: Obtain GID lazily, only when needed?
function rcall(f::Callable, p::Union{Integer,FunRef}; R::Type=eltype(f))
    rpar(R=R, f, p)[]::R
end

# macro rpar(args...)
#     if length(args) == 2
#         (p, expr) = args
#         resc(:(rpar(()->$expr, $p)))
#     elseif length(args) == 3
#         (R, p, expr) = args
#         resc(:(rpar(R=$R, ()->$expr, $p)))
#     else
#         error("Expected: @rpar [<type>] <proc> <expr>)")
#     end
# end

# macro rcall(args...)
#     if length(args) == 2
#         (p, expr) = args
#         resc(:(rcall(()->$expr, $p)))
#     elseif length(args) == 3
#         (R, p, expr) = args
#         resc(:(rcall(R=$R, ()->$expr, $p)))
#     else
#         error("Expected: @rcall [<type>] <proc> <expr>)")
#     end
# end

function make_local{R}(ref::FunRef{R})
    rpar(R=R, ref) do
        ref[]
    end::FunRef{R}
end
function get_local{R}(ref::FunRef{R})
    make_local(ref)[]::R
end

function unwrap(ref::FunRef{Any})
    r = FunRef{Any}()
    rgid = makegid(r)
    @schedule rexec(proc(ref)) do
        unwrap2(rgid, ref)
    end
    r::FunRef{Any}
end
function unwrap2(rgid::GID, ref::FunRef{Any})
    ref2 = ref[]::FunRef{Any}
    wait(ref2)
    rexec(rgid.proc) do
        unwrap3(rgid, ref2)
    end
end
function unwrap{R}(ref::FunRef{FunRef{R}})
    r = FunRef{R}()
    rgid = makegid(r)
    @schedule rexec(proc(ref)) do
        unwrap2(rgid, ref)
    end
    r::FunRef{R}
end
function unwrap2{R}(rgid::GID, ref::FunRef{FunRef{R}})
    ref2 = ref[]::FunRef{R}
    wait(ref2)
    rexec(rgid.proc) do
        unwrap3(rgid, ref2)
    end
end
function unwrap3{R}(rgid::GID, ref2::FunRef{R})
    r = getobj(rgid)::FunRef{R}
    free(rgid)
    set_from_ref!(r, ref2)
end

function remote(f::Callable, ref::FunRef; R::Type=eltype(f))
    r = FunRef{R}()
    rgid = makegid(r)
    @schedule rexec(proc(ref)) do
        remote2(R, rgid, f)
    end
    r::FunRef{R}
end
remote(f::Callable, p::Integer; R::Type=eltype(f)) = remote(f, Int(p), R=R)
function remote(f::Callable, p::Int; R::Type=eltype(f))
    r = FunRef{R}()
    rgid = makegid(r)
    rexec(p) do
        remote2(R, rgid, f)
    end
    r::FunRef{R}
end
function remote2(R::Type, rgid::GID, f::Callable)
    ref = FunRef{R}(f())
    rexec(rgid.proc) do
        remote3(rgid, Item(ref))
    end
end
function remote3(rgid::GID, ref::Item)
    r = getobj(rgid)
    free(rgid)
    set_from_ref!(r, ref.item::FunRef)
end

# macro remote(args...)
#     if length(args) == 2
#         (p, expr) = args
#         resc(:(remote(()->$expr, $p)))
#     elseif length(args) == 3
#         (R, p, expr) = args
#         resc(:(remote(R=$R, ()->$expr, $p)))
#     else
#         error("Expected: @remote [<type>] <proc> <expr>)")
#     end
# end

# FunRef is Foldable, Functor, Applicative, and Monad

@generated function freduce(op::Callable, zero, ref::FunRef, refs::FunRef...;
        R::Type=eltype(op))
    quote
        rcall(R=R, ref) do
            lrefs = ntuple(i->make_local(refs[i]), length(refs))
            op(zero, ref[], $([:(lrefs[$i][]) for i in 1:length(refs)]...))
        end
    end
end

# @generated function fmap(f::Callable, ref::FunRef, refs::FunRef...;
#         R::Type=eltype(f))
#     quote
#         remote(R=R, ref) do
#             lrefs = tuple($([:(make_local(refs[$i]))
#                              for i in 1:length(refs)]...))
#             f(ref[], $([:(lrefs[$i][]) for i in 1:length(refs)]...))
#         end
#     end
# end

function fmap(f::Callable, ref::FunRef, refs::FunRef...; R::Type=eltype(f))
    remote(R=R, ref) do
        lrefs = ntuple(i->make_local(refs[i]), length(refs))
        f(ref[], ntuple(i->lrefs[i][], length(refs))...)
    end
end

tycon{T,R}(::Type{FunRef{T}}, ::Type{R}) = FunRef{R}
valtype{T}(::Type{FunRef{T}}) = T

munit{T}(::Type{FunRef{T}}, x) = FunRef{T}(x)
mjoin{T}(xss::FunRef{FunRef{T}}) = unwrap(xss)
function mbind{T}(f::Callable, xs::FunRef{T}; R::Type=eltype(f))
    mjoin(fmap(R=R, f, xs))
end

end


text:
include("AIHelp.jl")
module BM
type Board

  state::Dict{Tuple{Int64,Int64},Tuple{String,Char}}
  currentPlayer::Char
  winner::Char # b, w, x for no winner

  Board(state, currentPlayer) = new(state, currentPlayer, 'x')

end
function startTestGame()
  board = Dict()
  board[(3,3)] = ("king",'b')
  board[(3,4)] = ("king",'w')
  return Board(board, 'b')
end

function startGame(gameType)
  board = Dict()
  if gameType == "standard"
    board[(5,1)] = ("king",'w')
    board[(5,9)] = ("king",'b')
    board[(4,1)] = ("gold general",'w')
    board[(6,1)] = ("gold general",'w')
    board[(4,9)] = ("gold general",'b')
    board[(6,9)] = ("gold general",'b')
    board[(3,1)] = ("silver general",'w')
    board[(7,1)] = ("silver general",'w')
    board[(3,9)] = ("silver general",'b')
    board[(7,9)] = ("silver general",'b')
    board[(2,1)] = ("knight",'w')
    board[(8,1)] = ("knight",'w')
    board[(2,9)] = ("knight",'b')
    board[(8,9)] = ("knight",'b')
    board[(1,1)] = ("lance",'w')
    board[(9,1)] = ("lance",'w')
    board[(1,9)] = ("lance",'b')
    board[(9,9)] = ("lance",'b')
    board[(2,2)] = ("bishop",'w')
    board[(8,8)] = ("bishop",'b')
    board[(8,2)] = ("rook",'w')
    board[(2,8)] = ("rook",'b')
    for i in 1:9
      board[(i,3)] = ("pawn",'w')
      board[(i,7)] = ("pawn",'b')
    end

  elseif gameType=="mini"

    board[(1,1)] = ("king",'w')
    board[(5,5)] = ("king",'b')
    board[(2,1)] = ("gold general",'w')
    board[(4,5)] = ("gold general",'b')
    board[(3,1)] = ("silver general",'w')
    board[(3,5)] = ("silver general",'b')
    board[(4,1)] = ("bishop",'w')
    board[(2,5)] = ("bishop",'b')
    board[(5,1)] = ("rook",'w')
    board[(1,5)] = ("rook",'b')
    board[(1,2)] = ("pawn",'w')
    board[(5,4)] = ("pawn",'b')

  elseif gameType=="chu"
    board[(1,1)] = ("lance",'w')
    board[(1,12)] = ("lance",'b')
    board[(12,1)] = ("lance",'w')
    board[(12,12)] = ("lance",'b')
    board[(2,1)] = ("ferocious leopard",'w')
    board[(2,12)] = ("ferocious leopard",'b')
    board[(11,1)] = ("ferocious leopard",'w')
    board[(11,12)] = ("ferocious leopard",'b')
    board[(3,1)] = ("copper general",'w')
    board[(3,12)] = ("copper general",'b')
    board[(10,1)] = ("copper general",'w')
    board[(10,12)] = ("copper general",'b')
    board[(4,1)] = ("silver general",'w')
    board[(4,12)] = ("silver general",'b')
    board[(9,1)] = ("silver general",'w')
    board[(9,12)] = ("silver general",'b')
    board[(5,1)] = ("gold general",'w')
    board[(5,12)] = ("gold general",'b')
    board[(8,1)] = ("gold general",'w')
    board[(8,12)] = ("gold general",'b')
    board[(6,1)] = ("king",'w')
    board[(7,12)] = ("king",'b')
    board[(7,1)] = ("drunk elephant",'w')
    board[(6,12)] = ("drunk elephant",'b')
    board[(1,2)] = ("reverse chariot",'w')
    board[(1,11)] = ("reverse chariot",'b')
    board[(12,2)] = ("reverse chariot",'w')
    board[(12,11)] = ("reverse chariot",'b')
    board[(3,2)] = ("bishop",'w')
    board[(3,11)] = ("bishop",'b')
    board[(10,2)] = ("bishop",'w')
    board[(10,11)] = ("bishop",'b')
    board[(5,2)] = ("blind tiger",'w')
    board[(5,11)] = ("blind tiger",'b')
    board[(8,2)] = ("blind tiger",'w')
    board[(8,11)] = ("blind tiger",'b')
    board[(6,2)] = ("kirin",'w')
    board[(7,11)] = ("kirin",'b')
    board[(7,2)] = ("phoenix",'w')
    board[(6,11)] = ("phoenix",'b')
    board[(1,3)] = ("side mover",'w')
    board[(1,10)] = ("side mover",'b')
    board[(12,3)] = ("side mover",'w')
    board[(12,10)] = ("side mover",'b')
    board[(2,3)] = ("vertical mover",'w')
    board[(2,10)] = ("vertical mover",'b')
    board[(11,3)] = ("vertical mover",'w')
    board[(11,10)] = ("vertical mover",'b')
    board[(3,3)] = ("rook",'w')
    board[(3,10)] = ("rook",'b')
    board[(10,3)] = ("rook",'w')
    board[(10,10)] = ("rook",'b')
    board[(4,3)] = ("dragon horse",'w')
    board[(4,10)] = ("dragon horse",'b')
    board[(9,3)] = ("dragon horse",'w')
    board[(9,10)] = ("dragon horse",'b')
    board[(5,3)] = ("dragon king",'w')
    board[(5,10)] = ("dragon king",'b')
    board[(8,3)] = ("dragon king",'w')
    board[(8,10)] = ("dragon king",'b')
    board[(6,3)] = ("lion",'w')
    board[(7,10)] = ("lion",'b')
    board[(7,3)] = ("queen",'w')
    board[(6,10)] = ("queen",'b')
    board[(4,5)] = ("go-between",'w')
    board[(4,8)] = ("go-between",'b')
    board[(9,5)] = ("go-between",'w')
    board[(9,8)] = ("go-between",'b')
    for j in 1:12
      board[(j,4)] = ("pawn",'w')
      board[(j,9)] = ("pawn",'b')
    end

  elseif gameType=="tenjiku"
    board[(1,1)] = ("lance",'w')
    board[(16,1)] = ("lance",'w')
    board[(2,1)] = ("knight",'w')
    board[(15,1)] = ("knight",'w')
    board[(3,1)] = ("ferocious leopard",'w')
    board[(14,1)] = ("ferocious leopard",'w')
    board[(4,1)] = ("iron general",'w')
    board[(13,1)] = ("iron general",'w')
    board[(5,1)] = ("copper general",'w')
    board[(12,1)] = ("copper general",'w')
    board[(6,1)] = ("silver general",'w')
    board[(11,1)] = ("silver general",'w')
    board[(7,1)] = ("gold general",'w')
    board[(10,1)] = ("gold general",'w')
    board[(9,1)] = ("drunk elephant",'w')
    board[(8,1)] = ("king",'w')
    board[(1,2)] = ("reverse chariot",'w')
    board[(16,2)] = ("reverse chariot",'w')
    board[(3,2)] = ("chariot soldier",'w')
    board[(14,2)] = ("chariot soldier",'w')
    board[(4,2)] = ("chariot soldier",'w')
    board[(13,2)] = ("chariot soldier",'w')
    board[(6,2)] = ("blind tiger",'w')
    board[(11,2)] = ("blind tiger",'w')
    board[(10,2)] = ("phoenix",'w')
    board[(7,2)] = ("kirin",'w')
    board[(9,2)] = ("queen",'w')
    board[(8,2)] = ("lion",'w')
    board[(1,3)] = ("side soldier",'w')
    board[(16,3)] = ("side soldier",'w')
    board[(2,3)] = ("vertical soldier",'w')
    board[(15,3)] = ("vertical soldier",'w')
    board[(3,3)] = ("bishop",'w')
    board[(14,3)] = ("bishop",'w')
    board[(4,3)] = ("dragon horse",'w')
    board[(13,3)] = ("dragon horse",'w')
    board[(5,3)] = ("dragon king",'w')
    board[(12,3)] = ("dragon king",'w')
    board[(6,3)] = ("water buffalo",'w')
    board[(11,3)] = ("water buffalo",'w')
    board[(7,3)] = ("fire demon",'w')
    board[(10,3)] = ("fire demon",'w')
    board[(9,3)] = ("free eagle",'w')
    board[(8,3)] = ("lion hawk",'w')
    board[(1,4)] = ("side mover",'w')
    board[(16,4)] = ("side mover",'w')
    board[(2,4)] = ("vertical mover",'w')
    board[(15,4)] = ("vertical mover",'w')
    board[(3,4)] = ("rook",'w')
    board[(14,4)] = ("rook",'w')
    board[(4,4)] = ("horned falcon",'w')
    board[(13,4)] = ("horned falcon",'w')
    board[(5,4)] = ("soaring eagle",'w')
    board[(12,4)] = ("soaring eagle",'w')
    board[(6,4)] = ("bishop general",'w')
    board[(11,4)] = ("bishop general",'w')
    board[(7,4)] = ("rook general",'w')
    board[(10,4)] = ("rook general",'w')
    board[(9,4)] = ("vice general",'w')
    board[(8,4)] = ("great general",'w')
    board[(5,6)] = ("dog",'w')
    board[(12,6)] = ("dog",'w')

    board[(1,16)] = ("lance",'b')
    board[(16,16)] = ("lance",'b')
    board[(2,16)] = ("knight",'b')
    board[(15,16)] = ("knight",'b')
    board[(3,16)] = ("ferocious leopard",'b')
    board[(14,16)] = ("ferocious leopard",'b')
    board[(4,16)] = ("iron general",'b')
    board[(13,16)] = ("iron general",'b')
    board[(5,16)] = ("copper general",'b')
    board[(12,16)] = ("copper general",'b')
    board[(6,16)] = ("silver general",'b')
    board[(11,16)] = ("silver general",'b')
    board[(7,16)] = ("gold general",'b')
    board[(10,16)] = ("gold general",'b')
    board[(8,16)] = ("drunk elephant",'b')
    board[(9,16)] = ("king",'b')
    board[(1,15)] = ("reverse chariot",'b')
    board[(16,15)] = ("reverse chariot",'b')
    board[(3,15)] = ("chariot soldier",'b')
    board[(14,15)] = ("chariot soldier",'b')
    board[(4,15)] = ("chariot soldier",'b')
    board[(13,15)] = ("chariot soldier",'b')
    board[(6,15)] = ("blind tiger",'b')
    board[(11,15)] = ("blind tiger",'b')
    board[(7,15)] = ("phoenix",'b')
    board[(10,15)] = ("kirin",'b')
    board[(8,15)] = ("queen",'b')
    board[(9,15)] = ("lion",'b')
    board[(1,14)] = ("side soldier",'b')
    board[(16,14)] = ("side soldier",'b')
    board[(2,14)] = ("vertical soldier",'b')
    board[(15,14)] = ("vertical soldier",'b')
    board[(3,14)] = ("bishop",'b')
    board[(14,14)] = ("bishop",'b')
    board[(4,14)] = ("dragon horse",'b')
    board[(13,14)] = ("dragon horse",'b')
    board[(5,14)] = ("dragon king",'b')
    board[(12,14)] = ("dragon king",'b')
    board[(6,14)] = ("water buffalo",'b')
    board[(11,14)] = ("water buffalo",'b')
    board[(7,14)] = ("fire demon",'b')
    board[(10,14)] = ("fire demon",'b')
    board[(8,14)] = ("free eagle",'b')
    board[(9,14)] = ("lion hawk",'b')
    board[(1,13)] = ("side mover",'b')
    board[(16,13)] = ("side mover",'b')
    board[(2,13)] = ("vertical mover",'b')
    board[(15,13)] = ("vertical mover",'b')
    board[(3,13)] = ("rook",'b')
    board[(14,13)] = ("rook",'b')
    board[(4,13)] = ("horned falcon",'b')
    board[(13,13)] = ("horned falcon",'b')
    board[(5,13)] = ("soaring eagle",'b')
    board[(12,13)] = ("soaring eagle",'b')
    board[(6,13)] = ("bishop general",'b')
    board[(11,13)] = ("bishop general",'b')
    board[(7,13)] = ("rook general",'b')
    board[(10,13)] = ("rook general",'b')
    board[(8,13)] = ("vice general",'b')
    board[(9,13)] = ("great general",'b')
    board[(12,11)] = ("dog",'b')
    board[(5,11)] = ("dog",'b')

    for k in 1:12
      board[(k,5)]=("pawn",'w')
      board[(k,12)]=("pawn",'b')
    end

  end
  return Board(board, 'b')
end

function currentPlayer(stateHist)
    return iseven(size(stateHist)[1]) ? 'w' : 'b' #size=1 is the starting state, no moves have been made
end

function nextState(board::Board, play::Tuple{Int64,Int64,Int64,Int64})
    newBoard = deepcopy(board)

    if ! ((play[3], play[4]) in newBoard.state.keys) #Target is empty
        newBoard.state[(play[3], play[4])] = newBoard.state[(play[1], play[2])]
        delete!(newBoard.state[(play[1], play[2])])
    else # target cell is not empty
        if (newBoard.state[(play[3], play[4])])[1] == "king"
            newBoard.winner = newBoard.currentPlayer
        end
        newBoard.state[(play[3], play[4])] = newBoard.state[(play[1], play[2])]
        delete!(newBoard,(play[1], play[2]))
    end

    newBoard.currentPlayer = (newBoard.currentPlayer == 'b' ? 'w' : 'b')

    return newBoard
end

function next_state(target, play::Tuple{Int64,Int64,Int64,Int64})
    newBoard = deepcopy(target)

    if ! ((play[3], play[4]) in newBoard.keys) #Target is empty
        newBoard[(play[3], play[4])] = newBoard[(play[1], play[2])]
        delete!(newBoard,(play[1], play[2]))
    else # target cell is not empty
        if (get(newBoard,(play[3], play[4]),"x"))[1] == "king"
            #newBoard.winner = newBoard.currentPlayer
        end
        newBoard[(play[3], play[4])] = newBoard[(play[1], play[2])]
        delete!(newBoard,(play[1], play[2]))
    end

    #newBoard.currentPlayer = (newBoard.currentPlayer == 'b' ? 'w' : 'b')

    return newBoard
end

function legalPlays(board::Board)
    return legalMovesPlayer(board.state, board.currentPlayer)
end

function winner(state)
  if !( ("king",'b') in values(state)) #=Hist[size(stateHist)[1]]) )=#
    return 'w'
  elseif !( ("king",'w') in values(state))
    return 'b'
  else
    return '?'
  end
end

export winner,legalPlays,next_state,nextState,currentPlayer,startTestGame,startGame
end

max_stars_repo_name: kinoroy/Shogi

max_issues_repo_name: kinoroy/Shulia

max_forks_repo_name: kinoroy/Shulia