text:
ys = cu([3. 3. 4. 4. 5.;
         5. 5. 6. 6. 7.])
us = cu(2*ones(2, 3, 4))
xs = CuArray{Int64}([1 2 3 4;
                     4 2 1 3;
                     3 5 5 3])

∇y_mul = [4. 4. 16. 4. 4.; 4. 4. 16. 4. 4.]
∇y_div = [.25 .25 .0625 .25 .25; .25 .25 .0625 .25 .25]
∇u_mean = cat([.5 .5 .25; .5 .5 .25], [.5 .5 .5; .5 .5 .5],
              [.25 .5 .5; .25 .5 .5], [.5 .25 .25; .5 .25 .25], dims=3)

@testset "cuda/grad" begin
    @testset "scatter" begin
        @test Zygote.gradient(x -> sum(scatter_add!(x, us, xs)), ys) == (ones(2, 5),)
        @test Zygote.gradient(x -> sum(scatter_add!(copy(ys), x, xs)), us) == (ones(2, 3, 4),)
        @test Zygote.gradient(x -> sum(scatter_add!(copy(ys), us, x)), xs) == (nothing,)

        @test Zygote.gradient(x -> sum(scatter_sub!(x, us, xs)), ys) == (ones(2, 5),)
        @test Zygote.gradient(x -> sum(scatter_sub!(copy(ys), x, xs)), us) == (-ones(2, 3, 4),)
        @test Zygote.gradient(x -> sum(scatter_sub!(copy(ys), us, x)), xs) == (nothing,)

        @test Zygote.gradient(x -> sum(scatter_max!(x, us, xs)), ys) == (ones(2, 5),)
        @test Zygote.gradient(x -> sum(scatter_max!(copy(ys), x, xs)), us) == (zeros(2, 3, 4),)
        @test Zygote.gradient(x -> sum(scatter_max!(copy(ys), us, x)), xs) == (nothing,)

        @test Zygote.gradient(x -> sum(scatter_min!(x, us, xs)), ys) == (zeros(2, 5),)
        @test Zygote.gradient(x -> sum(scatter_min!(copy(ys), x, xs)), us) == (ones(2, 3, 4),)
        @test Zygote.gradient(x -> sum(scatter_min!(copy(ys), us, x)), xs) == (nothing,)

        @test Zygote.gradient(x -> sum(scatter_mul!(x, us, xs)), ys) == (∇y_mul,)
        @test Zygote.gradient(x -> sum(scatter_mul!(copy(ys), x, xs)), us) == (2048*gather(ys, xs),)
        @test Zygote.gradient(x -> sum(scatter_mul!(copy(ys), us, x)), xs) == (nothing,)

        @test Zygote.gradient(x -> sum(scatter_div!(x, us, xs)), ys) == (∇y_div,)
        @test Zygote.gradient(x -> sum(scatter_div!(copy(ys), x, xs)), us) == (-gather(ys, xs)/8192,)
        @test Zygote.gradient(x -> sum(scatter_div!(copy(ys), us, x)), xs) == (nothing,)

        @test Zygote.gradient(x -> sum(scatter_mean!(x, us, xs)), ys) == (ones(2, 5),)
        @test Zygote.gradient(x -> sum(scatter_mean!(copy(ys), x, xs)), us) == (∇u_mean,)
        @test Zygote.gradient(x -> sum(scatter_mean!(copy(ys), us, x)), xs) == (nothing,)
    end
end


text:
@doc doc"""
    LinearAffineMetric <: Metric

The linear affine metric is the metric for symmetric positive definite matrices, that employs
matrix logarithms and exponentials, which yields a linear and affine metric.
"""
struct LinearAffineMetric <: RiemannianMetric end

is_default_metric(::SymmetricPositiveDefinite, ::LinearAffineMetric) = Val(true)

@doc doc"""
    distance(M::SymmetricPositiveDefinite, x, y)
    distance(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric})

Compute the distance on the [`SymmetricPositiveDefinite`](@ref) manifold between `x` and `y`,
as a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads

```math
d_{\mathcal P(n)}(x,y)
= \lVert \operatorname{Log}(x^{-\frac{1}{2}}yx^{-\frac{1}{2}})\rVert_{\mathrm{F}}.,
```
where $\operatorname{Log}$ denotes the matrix logarithm and
$\lVert\cdot\rVert_{\mathrm{F}}$ denotes the matrix Frobenius norm.
"""
function distance(M::SymmetricPositiveDefinite{N}, x, y) where {N}
    s = real.(eigvals(x, y))
    return any(s .<= eps()) ? 0 : sqrt(sum(abs.(log.(s)) .^ 2))
end

@doc doc"""
    exp(M::SymmetricPositiveDefinite, x, v)
    exp(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, x, v)

Compute the exponential map from `x` with tangent vector `v` on the
[`SymmetricPositiveDefinite`](@ref) `M` with its default [`MetricManifold`](@ref) having the
[`LinearAffineMetric`](@ref). The formula reads

```math
\exp_x v = x^{\frac{1}{2}}\operatorname{Exp}(x^{-\frac{1}{2}} v x^{-\frac{1}{2}})x^{\frac{1}{2}},
```

where $\operatorname{Exp}$ denotes to the matrix exponential.
"""
exp(::SymmetricPositiveDefinite, ::Any...)

function exp!(M::SymmetricPositiveDefinite{N}, y, x, v) where {N}
    e = eigen(Symmetric(x))
    U = e.vectors
    S = e.values
    Ssqrt = Diagonal(sqrt.(S))
    SsqrtInv = Diagonal(1 ./ sqrt.(S))
    xSqrt = Symmetric(U * Ssqrt * transpose(U))
    xSqrtInv = Symmetric(U * SsqrtInv * transpose(U))
    T = Symmetric(xSqrtInv * v * xSqrtInv)
    eig1 = eigen(T) # numerical stabilization
    Se = Diagonal(exp.(eig1.values))
    Ue = eig1.vectors
    xue = xSqrt * Ue
    return copyto!(y, xue * Se * transpose(xue))
end

@doc doc"""
    [Ξ,κ] = get_basis(M::SymmetricPositiveDefinite, x, B::DiagonalizingOrthonormalBasis)
    [Ξ,κ] = get_basis(M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric}, x, B::DiagonalizingOrthonormalBasis)

Return a orthonormal basis `Ξ` as a vector of tangent vectors (of length
[`manifold_dimension`](@ref) of `M`) in the tangent space of `x` on the
[`MetricManifold`](@ref) of [`SymmetricPositiveDefinite`](@ref) manifold `M` with
[`LinearAffineMetric`](@ref) that diagonalizes the curvature tensor $R(u,v)w$
with eigenvalues `κ` and where the direction `B.v` has curvature `0`.
"""
function get_basis(
    M::SymmetricPositiveDefinite{N},
    x,
    B::DiagonalizingOrthonormalBasis,
) where {N}
    xSqrt = sqrt(x)
    eigv = eigen(B.v)
    V = eigv.vectors
    Ξ = [
        (i == j ? 1 / 2 :
             1 / sqrt(2)) * (V[:, i] * transpose(V[:, j]) + V[:, j] * transpose(V[:, i]))
        for i = 1:N
        for j = i:N
    ]
    λ = eigv.values
    κ = [-1 / 4 * (λ[i] - λ[j])^2 for i = 1:N for j = i:N]
    return PrecomputedDiagonalizingOrthonormalBasis(Ξ, κ)
end
function get_basis(
    M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},
    x,
    B::DiagonalizingOrthonormalBasis,
) where {N}
    return get_basis(base_manifold(M), x, B)
end

function get_coordinates(
    M::SymmetricPositiveDefinite{N},
    x,
    v,
    B::ArbitraryOrthonormalBasis,
) where {N}
    dim = manifold_dimension(M)
    vout = similar(v, dim)
    @assert size(v) == (N, N)
    @assert dim == div(N * (N + 1), 2)
    k = 1
    for i = 1:N, j = i:N
        scale = ifelse(i == j, 1, sqrt(2))
        @inbounds vout[k] = v[i, j] * scale
        k += 1
    end
    return vout
end
function get_coordinates(
    M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},
    x,
    v,
    B::ArbitraryOrthonormalBasis,
) where {N}
    return get_coordinates(base_manifold(M), x, v, B)
end

function get_vector(
    M::SymmetricPositiveDefinite{N},
    x,
    v,
    B::ArbitraryOrthonormalBasis,
) where {N}
    dim = manifold_dimension(M)
    vout = allocate_result(M, get_vector, x)
    @assert size(v) == (div(N * (N + 1), 2),)
    @assert size(vout) == (N, N)
    k = 1
    for i = 1:N, j = i:N
        scale = ifelse(i == j, 1, 1 / sqrt(2))
        @inbounds vout[i, j] = v[k] * scale
        @inbounds vout[j, i] = v[k] * scale
        k += 1
    end
    return vout
end
function get_vector(
    M::MetricManifold{SymmetricPositiveDefinite{N},LinearAffineMetric},
    x,
    v,
    B::ArbitraryOrthonormalBasis,
) where {N}
    return get_vector(base_manifold(M), x, v, B)
end

@doc doc"""
    inner(M::SymmetricPositiveDefinite, x, v, w)
    inner(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, x, v, w)

Compute the inner product of `v`, `w` in the tangent space of `x` on
the [`SymmetricPositiveDefinite`](@ref) manifold `M`, as
a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads

````math
(v, w)_x = \operatorname{tr}(x^{-1} v x^{-1} w),
````
"""
function inner(M::SymmetricPositiveDefinite, x, v, w)
    F = cholesky(Symmetric(x))
    return tr((F \ Symmetric(v)) * (F \ Symmetric(w)))
end

@doc doc"""
    log(M::SymmetricPositiveDefinite, x, y)
    log(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, x, y)

Compute the logarithmic map from `x` to `y` on the [`SymmetricPositiveDefinite`](@ref)
as a [`MetricManifold`](@ref) with [`LinearAffineMetric`](@ref). The formula reads

```math
\log_x y =
x^{\frac{1}{2}}\operatorname{Log}(x^{-\frac{1}{2}}yx^{-\frac{1}{2}})x^{\frac{1}{2}},
```
where $\operatorname{Log}$ denotes to the matrix logarithm.
"""
log(::SymmetricPositiveDefinite, ::Any...)

function log!(M::SymmetricPositiveDefinite{N}, v, x, y) where {N}
    e = eigen(Symmetric(x))
    U = e.vectors
    S = e.values
    Ssqrt = Diagonal(sqrt.(S))
    SsqrtInv = Diagonal(1 ./ sqrt.(S))
    xSqrt = Symmetric(U * Ssqrt * transpose(U))
    xSqrtInv = Symmetric(U * SsqrtInv * transpose(U))
    T = Symmetric(xSqrtInv * y * xSqrtInv)
    e2 = eigen(T)
    Se = Diagonal(log.(max.(e2.values, eps())))
    xue = xSqrt * e2.vectors
    return mul!(v, xue, Se * transpose(xue))
end

@doc doc"""
    vector_transport_to(M::SymmetricPositiveDefinite, x, v, y, ::ParallelTransport)
    vector_transport_to(M::MetricManifold{SymmetricPositiveDefinite,LinearAffineMetric}, x, v, y, ::ParallelTransport)

Compute the parallel transport on the [`SymmetricPositiveDefinite`](@ref) as a
[`MetricManifold`](@ref) with the [`LinearAffineMetric`](@ref).
The formula reads

```math
P_{y\gets x}(v) = x^{\frac{1}{2}}
\operatorname{Exp}\bigl(
\frac{1}{2}x^{-\frac{1}{2}}\log_x(y)x^{-\frac{1}{2}}
\bigr)
x^{-\frac{1}{2}}v x^{-\frac{1}{2}}
\operatorname{Exp}\bigl(
\frac{1}{2}x^{-\frac{1}{2}}\log_x(y)x^{-\frac{1}{2}}
\bigr)
x^{\frac{1}{2}},
```

where $\operatorname{Exp}$ denotes the matrix exponential
and `log` the logarithmic map on [`SymmetricPositiveDefinite`](@ref)
(again with respect to the metric mentioned).
"""
vector_transport_to(::SymmetricPositiveDefinite, ::Any, ::Any, ::Any, ::ParallelTransport)

function vector_transport_to!(
    M::SymmetricPositiveDefinite{N},
    vto,
    x,
    v,
    y,
    ::ParallelTransport,
) where {N}
    distance(M, x, y) < 2 * eps(eltype(x)) && copyto!(vto, v)
    e = eigen(Symmetric(x))
    U = e.vectors
    S = e.values
    Ssqrt = sqrt.(S)
    SsqrtInv = Diagonal(1 ./ Ssqrt)
    Ssqrt = Diagonal(Ssqrt)
    xSqrt = Symmetric(U * Ssqrt * transpose(U))
    xSqrtInv = Symmetric(U * SsqrtInv * transpose(U))
    tv = Symmetric(xSqrtInv * v * xSqrtInv)
    ty = Symmetric(xSqrtInv * y * xSqrtInv)
    e2 = eigen(ty)
    Se = Diagonal(log.(e2.values))
    Ue = e2.vectors
    ty2 = Symmetric(Ue * Se * transpose(Ue))
    e3 = eigen(ty2)
    Sf = Diagonal(exp.(e3.values))
    Uf = e3.vectors
    xue = xSqrt * Uf * Sf * transpose(Uf)
    vtp = Symmetric(xue * tv * transpose(xue))
    return copyto!(vto, vtp)
end

ext: jl


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{typeof(Base.showarg),IOBuffer,OffsetArray{Int, 0, Array{Int, 0}},Bool})   # time: 0.037824474
    Base.precompile(Tuple{Type{IdOffsetRange{Int, Base.OneTo{Int}}},UnitRange{Int}})   # time: 0.009825722
    Base.precompile(Tuple{typeof(Base.inds2string),Tuple{IdOffsetRange{Int, Base.OneTo{Int}}, IdOffsetRange{Int, Base.OneTo{Int}}}})   # time: 0.0080779
    Base.precompile(Tuple{typeof(getindex),StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}},IdentityUnitRange{UnitRange{Int}}})   # time: 0.008055889
    Base.precompile(Tuple{typeof(getindex),StepRangeLen{Int, Int, Int},IdentityUnitRange{UnitRange{Int}}})   # time: 0.006584347
    Base.precompile(Tuple{typeof(getindex),StepRange{Int, Int},IdentityUnitRange{UnitRange{Int}}})   # time: 0.00650066
    Base.precompile(Tuple{typeof(getindex),LinRange{Float64},IdentityUnitRange{UnitRange{Int}}})   # time: 0.005844317
    Base.precompile(Tuple{typeof(zeros),Tuple{IdOffsetRange{Int, Base.OneTo{Int}}}})   # time: 0.007713056
    Base.precompile(Tuple{typeof(ones),Tuple{IdOffsetRange{Int, Base.OneTo{Int}}}})   # time: 0.007713056
    Base.precompile(Tuple{typeof(trues),Tuple{UnitRange{Int}, UnitRange{Int}}})   # time: 0.005478372
    Base.precompile(Tuple{typeof(falses),Tuple{UnitRange{Int}, UnitRange{Int}}})   # time: 0.005478372
    Base.precompile(Tuple{typeof(firstindex),IdOffsetRange{Int, Base.OneTo{Int}}})   # time: 0.004100289
end


text:
"""
    ForceIntensityModule

Module  to manage  distributed force intensity.
"""
module ForceIntensityModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.VectorCacheModule: VectorCache, updateretrieve!

"""
    ForceIntensity{T<:Number, F<:Function}

Distributed force (force intensity) type.
The force intensity class. The physical units are
force per unit volume, where volume depends on to which manifold
the force is applied:
- force/length^3 (when applied to a 3-D solid),
- force/length^2 (when applied to a surface),
- force/length^1 (when applied along a curve), or
- force/length^0 (when applied at a point).

Signature of the function to compute the value of the force  at any given point
`XYZ`, using the columns of the Jacobian matrix of the element, `tangents`,
and if necessary  also the finite element label, `fe_label`:

```
getforce!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
```

The buffer `forceout` is filled with the value  of the force.
"""
struct ForceIntensity{T<:Number, F<:Function}
    cache::VectorCache{T, F}   # cache of the current value of the force can be retrieved
end

"""
    ForceIntensity(::Type{T}, ndofn::FInt,
      computeforce!::F) where {T<:Number, F<:Function}

Construct force intensity when the function to compute the intensity vector is
given. This function needs to have a signature of
```
function computeforce!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
Calculate the force  and copy it into the buffer....
return forceout
end
```
and it needs to  fill in the buffer `forceout` with the current force at the
location `XYZ`, using if appropriate the information supplied in the Jacobian
matrix `tangents`, and the label of the finite element, `fe_label`.
"""
function ForceIntensity(::Type{T}, ndofn::FInt, computeforce!::F) where {T<:Number, F<:Function}
    # Allocate the buffer to be ready for the first call
    return ForceIntensity(VectorCache(computeforce!, zeros(T, ndofn)));
end

"""
    ForceIntensity(force::FVec{T}) where {T<:Number}

Construct force intensity when the constant force vector is given.
"""
function ForceIntensity(force::FVec{T}) where {T<:Number}
    function computeforce!(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        # do nothing:  the force is already in the buffer
        return forceout
    end
    return ForceIntensity(VectorCache(computeforce!, deepcopy(force)));
end

"""
    ForceIntensity(force::T) where {T<:Number}

Construct force intensity when the force is given as a scalar value.
The dimension of the force vector in this case is 1.
"""
function ForceIntensity(force::T) where {T<:Number}
    return ForceIntensity(T[force]);
end

"""
    updateforce!(self::ForceIntensity, ndofn::FInt, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)

Update the force intensity vector.

Returns a vector (stored in the cache `self.cache`).
"""
function updateforce!(self::ForceIntensity, XYZ::FFltMat,    tangents::FFltMat, fe_label::FInt)
    return updateretrieve!(self.cache, XYZ, tangents, fe_label)
end

end


text:
abstract type Computer    <: Machine
abstract type Workflow    <: Machine
abstract type Learner     <: Computer
abstract type Transformer <: Computer

function fit!(mc::Machine, input::DataFrame, output::Vector)
   error(typeof(mc)," has no implementation.")
end

function transform!(mc::Machine, input::DataFrame)
   error(typeof(mc)," has no implementation.")
end

function fit_transform!(mc::Machine, input::DataFrame, output::Vector)
   fit!(mc,input,output)
   transform!(mc,input)
end


text:
using SymbolicWedderburn
using PermutationGroups
using Cyclotomics

using DynamicPolynomials
using SumOfSquares
using SCS

C3 = PermGroup([perm"(1,2,3)"])
@polyvar x[1:3]

# @constraint m  sum(x.^2) - t in SOSCone() symmetry_group = C2
using MultivariateBases

mvec = monomials(x, 0:1)
dict = Dict(m => i for (i, m) in enumerate(mvec))

using MultivariatePolynomials
const MP = MultivariatePolynomials

function Base.:^(m::MP.AbstractMonomial, p::PermutationGroups.Perm)
    return prod(MP.variables(m).^(MP.exponents(m)^p))
end

G = PermGroup([Perm([dict[m] for m in mvec.^Ref(g)]) for g in gens(C3)])
ccG = conjugacy_classes(G)
chars = SymbolicWedderburn.characters_dixon(ccG)

c1, U1 = central_projection(chars[1])
c2, U2 = central_projection(chars[2])
c3, U3 = central_projection(chars[3])

@assert conj(U2) ==  U3

R1, ids1 = SymbolicWedderburn.row_echelon_form(float(U1))
R2, ids2 = SymbolicWedderburn.row_echelon_form(float(U2+U3))
R3, ids3 = SymbolicWedderburn.row_echelon_form(float(-E(4)*(U2-U3)))

msym = SOSModel(SCS.Optimizer)
@variable   msym t
@objective  msym Max t
@variable   msym sos1 SOSPoly(FixedPolynomialBasis(R1[1:length(ids1),:]*mvec))
@variable   msym sos2 SOSPoly(FixedPolynomialBasis(R2[1:length(ids2),:]*mvec))
@variable   msym sos3 SOSPoly(FixedPolynomialBasis(R3[1:length(ids3),:]*mvec))

@constraint msym sum(x.^2) - t == sos1 + sos2 + sos3
optimize!(msym)
#=
    Variables n = 10, constraints m = 19
    Cones:	primal zero / dual free vars: 10
	        soc vars: 9, soc blks: 3
=#

m = let
    m = SOSModel(SCS.Optimizer)
    @variable m t
    @objective m Max t
    @constraint m  sum(x.^2) - t in SOSCone()
    optimize!(m)

    #=
    Variables n = 11, constraints m = 20
    Cones:	primal zero / dual free vars: 10
	        sd vars: 10, sd blks: 1
    =#
end


text:
@testset "TrivialMStructure" begin
    b = StarAlgebras.Basis{UInt8}(words([:a, :b, :c], radius = 4))
    k = findfirst(w -> length(w) == 3, b) - 1

    mstr = StarAlgebras.TrivialMStructure{false}(b)

    @test mstr[1, 1] == 1
    @test mstr[1, 2] == 2
    @test mstr[2, 1] == 2
    @test mstr[3, 1] == 3

    idx = b[b[2]*b[3]] # == 8
    @test mstr[2, 3] == idx

    idx = b[b[3]*b[2]]
    @test mstr[3, 2] == idx

    tmstr = StarAlgebras.TrivialMStructure{true}(b)
    @test tmstr[1, 2] == 2
    @test tmstr[2, 3] == b[star(b[2])*b[3]]
    @test tmstr[3, 2] == b[star(b[3])*b[2]]

    @test sprint(show, MIME"text/plain"(), tmstr) ==
          "twisted TrivialMStructure over basis with $(length(basis(tmstr))) elements"

    @test_throws StarAlgebras.ProductNotDefined mstr[k+1, k]

    try
        w = mstr[k+1, k]
    catch ex
        @test ex isa StarAlgebras.ProductNotDefined
        @test sprint(Base.showerror, ex) ==
              "Product of elements 14 and 13 is not defined on the basis or the multiplicative structure could not be completed: a·a·a · c·c = a·a·a·c·c."
    end
end

@testset "MTable" begin
    b = StarAlgebras.Basis{UInt16}(words([:a, :b, :c, :d], radius = 4))
    k = findfirst(w -> length(w) == 3, b) - 1

    mstr = StarAlgebras.MTable(b, table_size = (k, k))

    @test mstr isa StarAlgebras.MTable{UInt16,false}

    @test all(mstr[i, i] ≠ 1 for i in 2:size(mstr, 1))
    @test all(mstr[1, i] == i for i in 1:size(mstr, 2))
    @test all(mstr[i, 1] == i for i in 1:size(mstr, 1))

    tmstr = StarAlgebras.MTable{true}(b, table_size = (k, k))

    @test tmstr isa StarAlgebras.MTable{UInt16,true}
    @test all(tmstr[i, i] != 1 for i in 2:size(tmstr, 1))
    @test all(tmstr[1, i] == i for i in 1:size(tmstr, 2))
    @test all(tmstr[i, 1] ≠ i for i in 1:size(tmstr, 1) if b[i] != star(b[i]))
end

@testset "CachedMTable" begin
    b = StarAlgebras.Basis{UInt8}(words([:a, :b, :c], radius = 4))
    k = findfirst(w -> length(w) == 3, b) - 1

    @test StarAlgebras.CachedMTable(b, table_size = (k, k)) isa
          StarAlgebras.CachedMTable{Word{Symbol},UInt8,typeof(b),Matrix{UInt8},false}

    @test StarAlgebras.CachedMTable{true}(b, table_size = (k, k)) isa
          StarAlgebras.CachedMTable{Word{Symbol},UInt8,typeof(b),Matrix{UInt8},true}

    @test StarAlgebras.CachedMTable(b, spzeros(UInt8, k, k)) isa StarAlgebras.CachedMTable{
        Word{Symbol},
        UInt8,
        typeof(b),
        SparseMatrixCSC{UInt8,Int64},
        false,
    }

    @test StarAlgebras.CachedMTable{true}(b, spzeros(UInt8, k, k)) isa
          StarAlgebras.CachedMTable{
        Word{Symbol},
        UInt8,
        typeof(b),
        SparseMatrixCSC{UInt8,Int64},
        true,
    }

    for mstr in [
        StarAlgebras.CachedMTable{false}(b, table_size = (k, k)),
        StarAlgebras.CachedMTable{true}(b, spzeros(UInt8, k, k)),
    ]

        @test all(iszero, mstr.table)
        StarAlgebras.cache!(mstr, 1, 2)
        @test mstr.table[1, 2] == 2
        @test mstr.table[1, 1] == 0

        idx = StarAlgebras._istwisted(mstr) ? b[star(b[2])*b[3]] : b[b[2]*b[3]]

        @test mstr.table[2, 3] == 0
        @test mstr[2, 3] == idx
        @test mstr.table[2, 3] == idx

        @test mstr.table[1, 3] == 0
        @test mstr.table[1, 4] == 0
        StarAlgebras.cache!(mstr, [1], [3, 4])
        @test mstr.table[1, 3] == 3
        @test mstr.table[1, 4] == 4

        @test_throws StarAlgebras.ProductNotDefined mstr[k+1, k]
    end

    tmstr = StarAlgebras.CachedMTable{true}(b, table_size = (k, k))

    @test all(iszero, tmstr.table)
    @test tmstr[1, 2] == 2
    @test tmstr[2, 3] == b[star(b[2])*b[3]]
    @test tmstr[3, 2] == b[star(b[3])*b[2]]
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: qldb

using AWS.Compat
using AWS.UUIDs
"""
    CancelJournalKinesisStream()

Ends a given Amazon QLDB journal stream. Before a stream can be canceled, its current status must be ACTIVE. You can't restart a stream after you cancel it. Canceled QLDB stream resources are subject to a 7-day retention period, so they are automatically deleted after this limit expires.

# Required Parameters
- `name`: The name of the ledger.
- `streamId`: The unique ID that QLDB assigns to each QLDB journal stream.

"""

cancel_journal_kinesis_stream(name, streamId; aws_config::AWSConfig=global_aws_config()) = qldb("DELETE", "/ledgers/$(name)/journal-kinesis-streams/$(streamId)"; aws_config=aws_config)
cancel_journal_kinesis_stream(name, streamId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("DELETE", "/ledgers/$(name)/journal-kinesis-streams/$(streamId)", args; aws_config=aws_config)

"""
    CreateLedger()

Creates a new ledger in your AWS account.

# Required Parameters
- `Name`: The name of the ledger that you want to create. The name must be unique among all of your ledgers in the current AWS Region. Naming constraints for ledger names are defined in Quotas in Amazon QLDB in the Amazon QLDB Developer Guide.
- `PermissionsMode`: The permissions mode to assign to the ledger that you want to create.

# Optional Parameters
- `DeletionProtection`: The flag that prevents a ledger from being deleted by any user. If not provided on ledger creation, this feature is enabled (true) by default. If deletion protection is enabled, you must first disable it before you can delete the ledger using the QLDB API or the AWS Command Line Interface (AWS CLI). You can disable it by calling the UpdateLedger operation to set the flag to false. The QLDB console disables deletion protection for you when you use it to delete a ledger.
- `Tags`: The key-value pairs to add as tags to the ledger that you want to create. Tag keys are case sensitive. Tag values are case sensitive and can be null.
"""

create_ledger(Name, PermissionsMode; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers", Dict{String, Any}("Name"=>Name, "PermissionsMode"=>PermissionsMode); aws_config=aws_config)
create_ledger(Name, PermissionsMode, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name, "PermissionsMode"=>PermissionsMode), args)); aws_config=aws_config)

"""
    DeleteLedger()

Deletes a ledger and all of its contents. This action is irreversible. If deletion protection is enabled, you must first disable it before you can delete the ledger using the QLDB API or the AWS Command Line Interface (AWS CLI). You can disable it by calling the UpdateLedger operation to set the flag to false. The QLDB console disables deletion protection for you when you use it to delete a ledger.

# Required Parameters
- `name`: The name of the ledger that you want to delete.

"""

delete_ledger(name; aws_config::AWSConfig=global_aws_config()) = qldb("DELETE", "/ledgers/$(name)"; aws_config=aws_config)
delete_ledger(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("DELETE", "/ledgers/$(name)", args; aws_config=aws_config)

"""
    DescribeJournalKinesisStream()

Returns detailed information about a given Amazon QLDB journal stream. The output includes the Amazon Resource Name (ARN), stream name, current status, creation time, and the parameters of your original stream creation request.

# Required Parameters
- `name`: The name of the ledger.
- `streamId`: The unique ID that QLDB assigns to each QLDB journal stream.

"""

describe_journal_kinesis_stream(name, streamId; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-kinesis-streams/$(streamId)"; aws_config=aws_config)
describe_journal_kinesis_stream(name, streamId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-kinesis-streams/$(streamId)", args; aws_config=aws_config)

"""
    DescribeJournalS3Export()

Returns information about a journal export job, including the ledger name, export ID, when it was created, current status, and its start and end time export parameters. This action does not return any expired export jobs. For more information, see Export Job Expiration in the Amazon QLDB Developer Guide. If the export job with the given ExportId doesn't exist, then throws ResourceNotFoundException. If the ledger with the given Name doesn't exist, then throws ResourceNotFoundException.

# Required Parameters
- `exportId`: The unique ID of the journal export job that you want to describe.
- `name`: The name of the ledger.

"""

describe_journal_s3_export(exportId, name; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-s3-exports/$(exportId)"; aws_config=aws_config)
describe_journal_s3_export(exportId, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-s3-exports/$(exportId)", args; aws_config=aws_config)

"""
    DescribeLedger()

Returns information about a ledger, including its state and when it was created.

# Required Parameters
- `name`: The name of the ledger that you want to describe.

"""

describe_ledger(name; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)"; aws_config=aws_config)
describe_ledger(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)", args; aws_config=aws_config)

"""
    ExportJournalToS3()

Exports journal contents within a date and time range from a ledger into a specified Amazon Simple Storage Service (Amazon S3) bucket. The data is written as files in Amazon Ion format. If the ledger with the given Name doesn't exist, then throws ResourceNotFoundException. If the ledger with the given Name is in CREATING status, then throws ResourcePreconditionNotMetException. You can initiate up to two concurrent journal export requests for each ledger. Beyond this limit, journal export requests throw LimitExceededException.

# Required Parameters
- `ExclusiveEndTime`: The exclusive end date and time for the range of journal contents that you want to export. The ExclusiveEndTime must be in ISO 8601 date and time format and in Universal Coordinated Time (UTC). For example: 2019-06-13T21:36:34Z  The ExclusiveEndTime must be less than or equal to the current UTC date and time.
- `InclusiveStartTime`: The inclusive start date and time for the range of journal contents that you want to export. The InclusiveStartTime must be in ISO 8601 date and time format and in Universal Coordinated Time (UTC). For example: 2019-06-13T21:36:34Z  The InclusiveStartTime must be before ExclusiveEndTime. If you provide an InclusiveStartTime that is before the ledger's CreationDateTime, Amazon QLDB defaults it to the ledger's CreationDateTime.
- `RoleArn`: The Amazon Resource Name (ARN) of the IAM role that grants QLDB permissions for a journal export job to do the following:   Write objects into your Amazon Simple Storage Service (Amazon S3) bucket.   (Optional) Use your customer master key (CMK) in AWS Key Management Service (AWS KMS) for server-side encryption of your exported data.
- `S3ExportConfiguration`: The configuration settings of the Amazon S3 bucket destination for your export request.
- `name`: The name of the ledger.

"""

export_journal_to_s3(ExclusiveEndTime, InclusiveStartTime, RoleArn, S3ExportConfiguration, name; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/journal-s3-exports", Dict{String, Any}("ExclusiveEndTime"=>ExclusiveEndTime, "InclusiveStartTime"=>InclusiveStartTime, "RoleArn"=>RoleArn, "S3ExportConfiguration"=>S3ExportConfiguration); aws_config=aws_config)
export_journal_to_s3(ExclusiveEndTime, InclusiveStartTime, RoleArn, S3ExportConfiguration, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/journal-s3-exports", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ExclusiveEndTime"=>ExclusiveEndTime, "InclusiveStartTime"=>InclusiveStartTime, "RoleArn"=>RoleArn, "S3ExportConfiguration"=>S3ExportConfiguration), args)); aws_config=aws_config)

"""
    GetBlock()

Returns a block object at a specified address in a journal. Also returns a proof of the specified block for verification if DigestTipAddress is provided. For information about the data contents in a block, see Journal contents in the Amazon QLDB Developer Guide. If the specified ledger doesn't exist or is in DELETING status, then throws ResourceNotFoundException. If the specified ledger is in CREATING status, then throws ResourcePreconditionNotMetException. If no block exists with the specified address, then throws InvalidParameterException.

# Required Parameters
- `BlockAddress`: The location of the block that you want to request. An address is an Amazon Ion structure that has two fields: strandId and sequenceNo. For example: {strandId:\"BlFTjlSXze9BIh1KOszcE3\",sequenceNo:14}
- `name`: The name of the ledger.

# Optional Parameters
- `DigestTipAddress`: The latest block location covered by the digest for which to request a proof. An address is an Amazon Ion structure that has two fields: strandId and sequenceNo. For example: {strandId:\"BlFTjlSXze9BIh1KOszcE3\",sequenceNo:49}
"""

get_block(BlockAddress, name; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/block", Dict{String, Any}("BlockAddress"=>BlockAddress); aws_config=aws_config)
get_block(BlockAddress, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/block", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("BlockAddress"=>BlockAddress), args)); aws_config=aws_config)

"""
    GetDigest()

Returns the digest of a ledger at the latest committed block in the journal. The response includes a 256-bit hash value and a block address.

# Required Parameters
- `name`: The name of the ledger.

"""

get_digest(name; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/digest"; aws_config=aws_config)
get_digest(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/digest", args; aws_config=aws_config)

"""
    GetRevision()

Returns a revision data object for a specified document ID and block address. Also returns a proof of the specified revision for verification if DigestTipAddress is provided.

# Required Parameters
- `BlockAddress`: The block location of the document revision to be verified. An address is an Amazon Ion structure that has two fields: strandId and sequenceNo. For example: {strandId:\"BlFTjlSXze9BIh1KOszcE3\",sequenceNo:14}
- `DocumentId`: The unique ID of the document to be verified.
- `name`: The name of the ledger.

# Optional Parameters
- `DigestTipAddress`: The latest block location covered by the digest for which to request a proof. An address is an Amazon Ion structure that has two fields: strandId and sequenceNo. For example: {strandId:\"BlFTjlSXze9BIh1KOszcE3\",sequenceNo:49}
"""

get_revision(BlockAddress, DocumentId, name; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/revision", Dict{String, Any}("BlockAddress"=>BlockAddress, "DocumentId"=>DocumentId); aws_config=aws_config)
get_revision(BlockAddress, DocumentId, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/revision", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("BlockAddress"=>BlockAddress, "DocumentId"=>DocumentId), args)); aws_config=aws_config)

"""
    ListJournalKinesisStreamsForLedger()

Returns an array of all Amazon QLDB journal stream descriptors for a given ledger. The output of each stream descriptor includes the same details that are returned by DescribeJournalKinesisStream. This action returns a maximum of MaxResults items. It is paginated so that you can retrieve all the items by calling ListJournalKinesisStreamsForLedger multiple times.

# Required Parameters
- `name`: The name of the ledger.

# Optional Parameters
- `max_results`: The maximum number of results to return in a single ListJournalKinesisStreamsForLedger request. (The actual number of results returned might be fewer.)
- `next_token`: A pagination token, indicating that you want to retrieve the next page of results. If you received a value for NextToken in the response from a previous ListJournalKinesisStreamsForLedger call, you should use that value as input here.
"""

list_journal_kinesis_streams_for_ledger(name; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-kinesis-streams"; aws_config=aws_config)
list_journal_kinesis_streams_for_ledger(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-kinesis-streams", args; aws_config=aws_config)

"""
    ListJournalS3Exports()

Returns an array of journal export job descriptions for all ledgers that are associated with the current AWS account and Region. This action returns a maximum of MaxResults items, and is paginated so that you can retrieve all the items by calling ListJournalS3Exports multiple times. This action does not return any expired export jobs. For more information, see Export Job Expiration in the Amazon QLDB Developer Guide.

# Optional Parameters
- `max_results`: The maximum number of results to return in a single ListJournalS3Exports request. (The actual number of results returned might be fewer.)
- `next_token`: A pagination token, indicating that you want to retrieve the next page of results. If you received a value for NextToken in the response from a previous ListJournalS3Exports call, then you should use that value as input here.
"""

list_journal_s3_exports(; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/journal-s3-exports"; aws_config=aws_config)
list_journal_s3_exports(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/journal-s3-exports", args; aws_config=aws_config)

"""
    ListJournalS3ExportsForLedger()

Returns an array of journal export job descriptions for a specified ledger. This action returns a maximum of MaxResults items, and is paginated so that you can retrieve all the items by calling ListJournalS3ExportsForLedger multiple times. This action does not return any expired export jobs. For more information, see Export Job Expiration in the Amazon QLDB Developer Guide.

# Required Parameters
- `name`: The name of the ledger.

# Optional Parameters
- `max_results`: The maximum number of results to return in a single ListJournalS3ExportsForLedger request. (The actual number of results returned might be fewer.)
- `next_token`: A pagination token, indicating that you want to retrieve the next page of results. If you received a value for NextToken in the response from a previous ListJournalS3ExportsForLedger call, then you should use that value as input here.
"""

list_journal_s3_exports_for_ledger(name; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-s3-exports"; aws_config=aws_config)
list_journal_s3_exports_for_ledger(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers/$(name)/journal-s3-exports", args; aws_config=aws_config)

"""
    ListLedgers()

Returns an array of ledger summaries that are associated with the current AWS account and Region. This action returns a maximum of 100 items and is paginated so that you can retrieve all the items by calling ListLedgers multiple times.

# Optional Parameters
- `max_results`: The maximum number of results to return in a single ListLedgers request. (The actual number of results returned might be fewer.)
- `next_token`: A pagination token, indicating that you want to retrieve the next page of results. If you received a value for NextToken in the response from a previous ListLedgers call, then you should use that value as input here.
"""

list_ledgers(; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers"; aws_config=aws_config)
list_ledgers(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/ledgers", args; aws_config=aws_config)

"""
    ListTagsForResource()

Returns all tags for a specified Amazon QLDB resource.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) for which you want to list the tags. For example:  arn:aws:qldb:us-east-1:123456789012:ledger/exampleLedger

"""

list_tags_for_resource(resourceArn; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("GET", "/tags/$(resourceArn)", args; aws_config=aws_config)

"""
    StreamJournalToKinesis()

Creates a journal stream for a given Amazon QLDB ledger. The stream captures every document revision that is committed to the ledger's journal and delivers the data to a specified Amazon Kinesis Data Streams resource.

# Required Parameters
- `InclusiveStartTime`: The inclusive start date and time from which to start streaming journal data. This parameter must be in ISO 8601 date and time format and in Universal Coordinated Time (UTC). For example: 2019-06-13T21:36:34Z  The InclusiveStartTime cannot be in the future and must be before ExclusiveEndTime. If you provide an InclusiveStartTime that is before the ledger's CreationDateTime, QLDB effectively defaults it to the ledger's CreationDateTime.
- `KinesisConfiguration`: The configuration settings of the Kinesis Data Streams destination for your stream request.
- `RoleArn`: The Amazon Resource Name (ARN) of the IAM role that grants QLDB permissions for a journal stream to write data records to a Kinesis Data Streams resource.
- `StreamName`: The name that you want to assign to the QLDB journal stream. User-defined names can help identify and indicate the purpose of a stream. Your stream name must be unique among other active streams for a given ledger. Stream names have the same naming constraints as ledger names, as defined in Quotas in Amazon QLDB in the Amazon QLDB Developer Guide.
- `name`: The name of the ledger.

# Optional Parameters
- `ExclusiveEndTime`: The exclusive date and time that specifies when the stream ends. If you don't define this parameter, the stream runs indefinitely until you cancel it. The ExclusiveEndTime must be in ISO 8601 date and time format and in Universal Coordinated Time (UTC). For example: 2019-06-13T21:36:34Z
- `Tags`: The key-value pairs to add as tags to the stream that you want to create. Tag keys are case sensitive. Tag values are case sensitive and can be null.
"""

stream_journal_to_kinesis(InclusiveStartTime, KinesisConfiguration, RoleArn, StreamName, name; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/journal-kinesis-streams", Dict{String, Any}("InclusiveStartTime"=>InclusiveStartTime, "KinesisConfiguration"=>KinesisConfiguration, "RoleArn"=>RoleArn, "StreamName"=>StreamName); aws_config=aws_config)
stream_journal_to_kinesis(InclusiveStartTime, KinesisConfiguration, RoleArn, StreamName, name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/ledgers/$(name)/journal-kinesis-streams", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("InclusiveStartTime"=>InclusiveStartTime, "KinesisConfiguration"=>KinesisConfiguration, "RoleArn"=>RoleArn, "StreamName"=>StreamName), args)); aws_config=aws_config)

"""
    TagResource()

Adds one or more tags to a specified Amazon QLDB resource. A resource can have up to 50 tags. If you try to create more than 50 tags for a resource, your request fails and returns an error.

# Required Parameters
- `Tags`: The key-value pairs to add as tags to the specified QLDB resource. Tag keys are case sensitive. If you specify a key that already exists for the resource, your request fails and returns an error. Tag values are case sensitive and can be null.
- `resourceArn`: The Amazon Resource Name (ARN) to which you want to add the tags. For example:  arn:aws:qldb:us-east-1:123456789012:ledger/exampleLedger

"""

tag_resource(Tags, resourceArn; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/tags/$(resourceArn)", Dict{String, Any}("Tags"=>Tags); aws_config=aws_config)
tag_resource(Tags, resourceArn, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Tags"=>Tags), args)); aws_config=aws_config)

"""
    UntagResource()

Removes one or more tags from a specified Amazon QLDB resource. You can specify up to 50 tag keys to remove.

# Required Parameters
- `resourceArn`: The Amazon Resource Name (ARN) from which you want to remove the tags. For example:  arn:aws:qldb:us-east-1:123456789012:ledger/exampleLedger
- `tagKeys`: The list of tag keys that you want to remove.

"""

untag_resource(resourceArn, tagKeys; aws_config::AWSConfig=global_aws_config()) = qldb("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), args)); aws_config=aws_config)

"""
    UpdateLedger()

Updates properties on a ledger.

# Required Parameters
- `name`: The name of the ledger.

# Optional Parameters
- `DeletionProtection`: The flag that prevents a ledger from being deleted by any user. If not provided on ledger creation, this feature is enabled (true) by default. If deletion protection is enabled, you must first disable it before you can delete the ledger using the QLDB API or the AWS Command Line Interface (AWS CLI). You can disable it by calling the UpdateLedger operation to set the flag to false. The QLDB console disables deletion protection for you when you use it to delete a ledger.
"""

update_ledger(name; aws_config::AWSConfig=global_aws_config()) = qldb("PATCH", "/ledgers/$(name)"; aws_config=aws_config)
update_ledger(name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = qldb("PATCH", "/ledgers/$(name)", args; aws_config=aws_config)


text:
using NCDatasets
using Test
using Dates
using Printf

using Compat

sz = (4,5)
filename = tempname()
# The mode "c" stands for creating a new file (clobber)
ds = Dataset(filename,"c")

# define the dimension "lon" and "lat"
defDim(ds,"lon",sz[1])
defDim(ds,"lat",sz[2])

# variables
for T in [UInt8,Int8,UInt16,Int16,UInt32,Int32,UInt64,Int64,Float32,Float64]
    #for T in [Float32]
    local v

    # write array
    v = defVar(ds,"var-$T",T,("lon","lat"))

    # write array
    v[:,:] = fill(T(123),size(v))
    @test all(v[:,:][:] .== 123)

    # write scalar
    v[:,:] = T(100)
    @test all(v[:,:][:] .== 100)

    # write array (different type)
    v[:,:] = fill(123,size(v))
    @test all(v[:,:][:] .== 123)

    # write scalar (different type)
    v[:,:] = 100
    @test all(v[:,:][:] .== 100)

    # using StepRange as index
    # write array
    v[1:end,1:end] = fill(T(123),size(v))
    @test all(v[:,:][:] .== 123)

    # write scalar
    v[1:end,1:end] = T(100)
    @test all(v[:,:][:] .== 100)

    # write array (different type)
    v[1:end,1:end] = fill(123,size(v))
    @test all(v[:,:][:] .== 123)

    # write scalar (different type)
    v[1:end,1:end] = 100
    @test all(v[:,:][:] .== 100)

    # step range
    ref = zeros(sz)
    v[:,:] = 0

    ref[1:2:end,2:2:end] .= 1
    v[1:2:end,2:2:end] = 1
    @test v[:,:] == ref

    # write scalar (different type)
    ref[1:2:end,2:2:end] .= UInt8(2)
    v[1:2:end,2:2:end] = UInt8(2)
    @test v[:,:] == ref

    ref[1,1] = 3
    v[1,1] = 3

    @test v[:,:] == ref
end

v = defVar(ds,"var-Char",Char,("lon","lat"))

# write array (without transformation)
v.var[:,:] = fill('a',size(v))
@test all(i -> i == 'a',v.var[:,:][:])

v.var[1,1] = 'x'
@test v.var[1,1] == 'x'
@test v.var[1,2] == 'a'

# write scalar
v.var[:,:] .= 'b'
@test all(i -> i == 'b',v.var[:,:][:])

# write array (with transformation)
v[:,:] = fill('c',size(v))
@test all(i -> i == 'c',v.var[:,:][:])

# write scalar
v[:,:] = 'd'
@test all(i -> i == 'd',v.var[:,:][:])

# using StepRange as index
# write array (without transformation)
v.var[1:end,1:end] = fill('e',size(v))
@test all(i -> i == 'e',v.var[:,:][:])

# write scalar
v.var[1:end,1:end] = 'f'
@test all(i -> i == 'f',v.var[:,:][:])

# write array (with transformation)
v[1:end,1:end] = fill('g',size(v))
@test all(i -> i == 'g',v.var[:,:][:])

# write scalar
v[1:end,1:end] = 'h'
@test all(i -> i == 'h',v.var[:,:][:])

# write with StepRange
v[:,:] = 'h'
ref = fill('h',sz)

ref[1:2:end,1:2:end] .= Ref('i')
#ref[1:2:end,1:2:end] .= 'i'
v[1:2:end,1:2:end] = 'i'

@test v[:,:] == ref

ref = ['a'+i+2*j for i = 0:sz[1]-1, j= 0:sz[2]-1]
v[:,:] = ref

@test v.var[1:1,1:1] == ref[1:1,1:1]
@test v.var[1,1] == ref[1,1]
@test v.var[2:3,1] == ref[2:3,1]

#@test v[1,1] == ref[1,1]

@test dimnames(v) == ("lon","lat")
@test name(v) == "var-Char"

close(ds)


text:
#=
# 4. Multiple stationary bodies
Adding multiple bodies to a problem is easy, using the concepts of a `BodyList`
and `RigidTransformList`.
=#
using ViscousFlow
#-
using Plots

#=
In this example, we will set up a problem with three cylinders arranged in a
formation in a free stream.
=#
Re = 200 # Reynolds number
U = 1.0 # Free stream velocity
U∞ = (U,0.0);

#=
Set up the domain, grid spacing, and time step size
=#
xlim = (-2.0,4.0)
ylim = (-2.0,2.0)
Δx, Δt = setstepsizes(Re,gridRe=4)

#=
### Set up bodies
We start by initializing a `BodyList` and an associated `RigidTransformList`.
Each member of the `RigidTransformList` will be used to place the respective body in
the correct position and orientation.
=#
bl = BodyList()
tl = RigidTransformList()

# Place the first body at (-1,0)
push!(bl,Circle(0.5,1.5Δx))
push!(tl,RigidTransform((-1.,0.),0.));

# Place the second body at (1,-1)
push!(bl,Circle(0.5,1.5Δx))
push!(tl,RigidTransform((1.,-1.),0.));

# and place the third body at (1,1)
push!(bl,Circle(0.5,1.5Δx))
push!(tl,RigidTransform((1.,1.),0.))

# Perform the actual transformation. Note that this operation works `in-place`:
tl(bl)

# #### Plot the initial configuration of the bodies
# Just to check they are in the right places
plot(bl,xlim=xlim,ylim=ylim)

#=
### Construct the system structure
We construct the system with the same syntax as for a single body:
=#
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,bl,freestream = U∞)
#-
u0 = newstate(sys)
tspan = (0.0,10.0)
integrator = init(u0,tspan,sys)
#=
### Solve
Here, we run it for only a little while, just to demonstrate:
=#
@time step!(integrator,0.5)
#=
### Examine
Let's make an animation
=#
sol = integrator.sol;
@gif for (u,t) in zip(sol.u,sol.t)
    plot(vorticity(u,sys,t),sys,clim=(-10,10),levels=range(-10,10,length=30), color = :RdBu)
end every 5

# Now we will examine the force on each body
fx1, fy1 = force(sol,sys,1)
fx2, fy2 = force(sol,sys,2)
fx3, fy3 = force(sol,sys,3);
#-
plt = plot(layout = (2,1), size = (600, 400))
plot!(plt[1],sol.t,2*fx1,xlim=(0,20),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",label="Lead body",title="Drag force")
plot!(plt[2],sol.t,2*fy1,xlim=(0,20),ylim=(-2,2),xlabel="Convective time",ylabel="\$C_L\$",label="Lead body",title="Side force")
plot!(plt[1],sol.t,2*fx2,xlim=(0,20),ylim=(0,4),xlabel="Convective time",ylabel="\$C_D\$",label="Trailing body",title="Drag force")
plot!(plt[2],sol.t,2*fy2,xlim=(0,20),ylim=(-2,2),xlabel="Convective time",ylabel="\$C_L\$",label="Trailing body",title="Side force")
#-
println("Mean drag coefficient on lead body = ", GridUtilities.mean(2*fx1))
#-
println("Mean drag coefficient on trailing body = ", GridUtilities.mean(2*fx2))


text:
# This tutorial is loosely adapted from [this pandas tutorial](https://pandas.pydata.org/pandas-docs/stable/getting_started/10min.html) as well as the [DataFrames.jl documentation](http://juliadata.github.io/DataFrames.jl/latest/man/getting_started/).
# It is by no means meant to be a complete introduction, rather, it focuses on some key functionalities that are particularly useful in a classical machine learning context.
#
# ## Basics
#
# To start with, we will use the `Boston` dataset which is very simple.

using RDatasets, DataFrames

boston = dataset("MASS", "Boston");

# The `dataset` function returns a `DataFrame` object:

typeof(boston)

# ### Accessing data
#
# Intuitively a DataFrame is just a wrapper around a number of columns, each of which is a `Vector` of some type with a name"

names(boston)

# You can view the first few rows using `first` and specifying a number of rows:

first(boston, 4)

# You can access one of those columns easily using `.colname`, this returns a vector that you can access like any Julia vector:

boston.Crim[1:5]

# You can also just access the dataframe as you would a big matrix:

boston[3, 5]

# or specifying a range of rows/columns:

boston[1:5, [:Crim, :Zn]]

# or, similarly,

boston[1:5, 1:2]

# The `select` function is very convenient to get sub dataframes of interest:

b1 = select(boston, [:Crim, :Zn, :Indus])
first(b1, 2)

# The `Not` syntax is  also very  useful:

b2 = select(boston, Not(:NOx))
first(b2, 2)

# Finally, if you would like to drop columns, you can use `select!` which will mutate the dataframe in place:

select!(b1, Not(:Crim))
first(b1, 2)

# ### Describing the data
#
# `StatsBase` offers a convenient `describe` function which you can use on a DataFrame to get an overview of the data:

using StatsBase
describe(boston, :min, :max, :mean, :median, :std)

# You can pass a number of symbols to the `describe` function to indicate which statistics to compute for each feature:
#
# * `mean`, `std`, `min`, `max`, `median`, `first`, `last` are all fairly self explanatory
# * `q25`, `q75` are respectively for the 25th and 75th percentile,
# * `eltype`, `nunique`, `nmissing` can also be used
#
# You can also  pass your custom function with a pair `name => function` for instance:

foo(x) = sum(abs.(x)) / length(x)
d = describe(boston, :mean, :median, :foo => foo)
first(d, 3)

# The `describe` function returns a derived object with one row per feature and one column per required statistic.
#
# Further to `StatsBase`, `Statistics` offers a range of useful functions for data analysis.

using Statistics

#
# ### Converting the data
#
# If you want to get the content of the dataframe as one big matrix, use `convert`:

mat = convert(Matrix, boston)
mat[1:3, 1:3]

# ### Adding columns
#
# Adding a column to a dataframe is very easy:

boston.Crim_x_Zn = boston.Crim .* boston.Zn;

# that's it! Remember also that you can drop columns or make subselections with `select` and `select!`.

# ### Missing values
#
# Let's load a dataset with missing values

mao = dataset("gap", "mao")
describe(mao, :nmissing)

# Lots of missing values...
# If  you wanted to compute simple functions on columns, they  may just return `missing`:

std(mao.Age)

# The `skipmissing` function can help counter this  easily:

std(skipmissing(mao.Age))

#
# ## Split-Apply-Combine
#
# This is a shorter version of the [DataFrames.jl tutorial](http://juliadata.github.io/DataFrames.jl/latest/man/split_apply_combine/).

iris = dataset("datasets", "iris")
first(iris, 3)

# ### `groupby`
#
# The `groupby` function allows to form "sub-dataframes" corresponding to groups of rows.
# This can be very convenient to run specific analyses for specific groups without copying the data.
#
# The basic usage is `groupby(df, cols)` where `cols` specifies one or several columns to use for the grouping.
#
# Consider a simple example: in `iris` there is a `Species` column with 3 species:

unique(iris.Species)

# We can form views for each of these:

gdf = groupby(iris, :Species);

# The `gdf` object now corresponds to **views** of the original dataframe for each of the 3 species; the first species is `"setosa"` with:

subdf_setosa = gdf[1]
describe(subdf_setosa, :min, :mean, :max)

# Note that `subdf_setosa` is a `SubDataFrame` meaning that it is just a view of the parent dataframe `iris`; if you modify that parent dataframe then the sub dataframe is also  modified.
#
# See `?groupby` for more information.

# ### `by`
#
# The `by` function allows you to compute statistics on given columns for each group.
# Let's start with a very simple example:

by(iris, :Species, :PetalLength => mean)

# So this returns the mean of the `:PetaLength` feature for each `:Species`.
#
# You can do this for several columns/statistics at the time and give specific names to the results:

by(iris, :Species, MPL = :PetalLength => mean, SPL = :PetalLength => std)

#
# See `?by` for more information.
#
# ### `aggregate`
#
# The aggregate function is a bit like `by` except that it applies the given function on all columns:

DataFrames.aggregate(iris, :Species, std)

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
# Generic instrument struct
abstract type Instrument end

mutable struct INSTR{instr_name} <: Instrument
    name::Symbol
    address::String
    handle# this needs to be anything, to accept PXI instruments
	initialized::Bool
	bufSize::UInt32
end
# Generic instrument constructor
INSTR(instr_name, address) = INSTR{instr_name}(instr_name, address, 0, false, UInt32(1024))

##############################################
# PSU
E364x1 = Union{INSTR{:KeysightE3645A},INSTR{:AgilentE3645A}}
E364x2 = Union{INSTR{:KeysightE3646A},INSTR{:AgilentE3646A},INSTR{:KeysightE3649A},INSTR{:AgilentE3649A}}

##############################################
# AWG
F332x1 = Union{INSTR{:Agilent33250A}}
F335x2 = Union{INSTR{:Agilent33500B},INSTR{:Keysight33500B},INSTR{:Agilent33522A}}
NIFgen = Union{INSTR{:NI5421},INSTR{:NI5422}}
DG1000Z = Union{INSTR{:RigolDG1022Z}}

##############################################
# DMM
KE200x = Union{INSTR{:Keithley2000},INSTR{:KE2000},INSTR{:Keithley2001},INSTR{:KE2001},}
NIDmm = Union{INSTR{:NI4071}}

##############################################
# Lecroy SCOPE
LCR6xxx = Union{INSTR{:HDO6054A},INSTR{:HRO64zi}}
# NI PXI SCOPE
NIScope = Union{INSTR{:NI5122},INSTR{:NI5114},INSTR{:NI5152},}