text:
"""
    IncRel{LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}

Incidence relation expressing connectivity between an entity of the shape on the
left and a list of entities of the shape on the right.

The incidence relation may be referred to by its code: `(d1, d2)`, where `d1` is
the manifold dimension of the shape collection on the left, and `d2` is the
manifold dimension of the shape collection on the right.

Main operations referred to with the abbreviations used in the paper:

- skt = ir_skeleton
- trp = ir_transpose
- bbf = ir_bbyfacets
- bbr = ir_bbyridges
- idt = ir_identity

Access to the incidence relation data is through indexing:
For instance
```
ir[55]
```
provides access to the entire incidence relation `ir` of the entity 55.
The first component of this incidence relation for the entity 55 is
```
ir[55][1]
```
or equivalently
```
ir[55, 1]
```
"""
struct IncRel{LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}
    left::ShapeColl{LEFT}
    right::ShapeColl{RIGHT}
    _v::Vector{T} # vector of vectors of shape numbers
    name::String # name of the incidence relation
    function IncRel(
        left::ShapeColl{LEFT},
        right::ShapeColl{RIGHT},
        v::Vector{T},
        name::String,
        internal
    ) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}
        @_check (nshapes(left) == length(v))
        minn = typemax(eltype(v[1]))
        maxn = typemin(eltype(v[1]))
        for i in 1:length(v)
            minn = min(minn, minimum(abs.(v[i])))
            maxn = max(maxn, maximum(abs.(v[i])))
        end
        @_check (nshapes(right) >= maxn)
        return new{LEFT,RIGHT,T}(left, right, v, name)
    end
end

# Provide access to the incidence relation data as if it was a one-dimensional
# or two dimensional array.
Base.IndexStyle(::Type{<:IncRel}) = IndexLinear()
Base.size(ir::IR) where {IR<:IncRel} = (length(ir._v),)
Base.getindex(ir::IR, i::Int) where {IR<:IncRel} = ir._v[i]
Base.getindex(ir::IR, i::Int, k::Int) where {IR<:IncRel} = ir._v[i][k]

"""
    IncRel(left::ShapeColl{LEFT}, right::ShapeColl{RIGHT}, v::Vector{T}) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}

Convenience constructor with a vector of vectors and a default name.
"""
function IncRel(
    left::ShapeColl{LEFT},
    right::ShapeColl{RIGHT},
    v::Vector{T}
) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}
    IncRel(left, right, deepcopy(v), "(" * left.name * ", " * right.name * ")", true)
end

"""
    IncRel(left::ShapeColl{LEFT}, right::ShapeColl{RIGHT}, v::Vector{T}, name::String) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}

Convenience constructor with a vector of vectors and a name.
"""
function IncRel(
    left::ShapeColl{LEFT},
    right::ShapeColl{RIGHT},
    v::Vector{T},
    name::String
) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}
    IncRel(left, right, deepcopy(v), name, true)
end

"""
    IncRel(left::ShapeColl{LEFT}, right::ShapeColl{RIGHT}, data::Matrix{MT}) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,MT}

Convenience constructor supplying a matrix instead of a vector of vectors and a default name.
"""
function IncRel(
    left::ShapeColl{LEFT},
    right::ShapeColl{RIGHT},
    data::Matrix{MT}
) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,MT}
    _v = [SVector{nvertices(left.shapedesc)}(data[idx, :]) for idx in 1:size(data, 1)]
    IncRel(left, right, _v, "(" * left.name * ", " * right.name * ")", true)
end

"""
    IncRel(left::ShapeColl{LEFT}, right::ShapeColl{RIGHT}, data::Matrix{MT}, name::String) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,MT}

Convenience constructor supplying a matrix instead of a vector of vectors and a name.
"""
function IncRel(
    left::ShapeColl{LEFT},
    right::ShapeColl{RIGHT},
    data::Matrix{MT},
    name::String
) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,MT}
    _v = [SVector{nvertices(left.shapedesc)}(data[idx, :]) for idx in 1:size(data, 1)]
    IncRel(left, right, _v, name, true)
end

"""
    indextype(ir::IncRel)

Provide type of the incidence index (entity number).
Presumably some type of integer.
"""
indextype(ir::IncRel) = eltype(eltype(ir._v))

"""
    nrelations(ir::IncRel)

Number of individual relations in the incidence relation.
"""
nrelations(ir::IncRel) = length(ir._v)

"""
    nentities(ir::IncRel, j::IT) where {IT}

Number of individual entities in the `j`-th relation in the incidence relation.
"""
nentities(ir::IncRel, j::IT) where {IT} = length(ir._v[j])

"""
    code(ir::IncRel{LEFT,RIGHT,T}) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T}

Formulate the code of the incidence relation.
"""
ir_code(ir::IncRel{LEFT,RIGHT,T}) where {LEFT<:AbsShapeDesc,RIGHT<:AbsShapeDesc,T} =
    (manifdim(shapedesc(ir.left)), manifdim(shapedesc(ir.right)))

"""
    ir_transpose(ir::IncRel, name = "trp")

Compute the "transposed" incidence relation `(d1, d2)`, where `d2 >= d1`.

This only makes sense for `d2 >= 0`. For `d2=1` we get for each vertex the list
of edges connected to the vertex, and analogously faces and cells for `d=2` and
`d=3`.

# Returns
Incidence relation for the transposed incidence relation. The left and right
shape collection are swapped in the output relative to the input.
"""
function ir_transpose(ir::IncRel, name = "trp")
    @_check (manifdim(ir.left) >= manifdim(ir.right))
    inttype = eltype(ir._v[1])
    # Find out how many of the transpose incidence relations there are
    nvmax = 0
    for j in 1:nrelations(ir)
        for k in 1:nentities(ir, j)
            nvmax = max(nvmax, ir[j, k])
        end
    end
    # pre-allocate relations vector
    _v = Vector{inttype}[]
    sizehint!(_v, nvmax)
    # Initialize the relations to empty
    for i in 1:nvmax
        push!(_v, inttype[])  # initially empty arrays
    end
    # Build the transpose relations
    for j in 1:nrelations(ir)
        for k in 1:nentities(ir, j)
            c = abs(ir[j, k]) # this could be an oriented entity: remove the sign
            push!(_v[c], j)
        end
    end
    return IncRel(ir.right, ir.left, _v, name)
end

function _asmatrix(ir, inttype)
    c = fill(inttype(0), nshapes(ir.left), nvertices(shapedesc(ir.left)))
    for i in 1:nshapes(ir.left)
        c[i, :] = ir[i]
    end
    return c
end

function _mysortdim2!(A)
    # Sort each row  of A in ascending order.
    m, n = size(A)
    r = zeros(eltype(A), n)
    for k in 1:m
        for i in 1:n
            r[i] = A[k, i]
        end
        sort!(r)
        for i in 1:n
            A[k, i] = r[i]
        end
    end
    return A
end

function _mysortrowsperm(A)
    # Sort the rows of A by sorting each column from back to front.
    m, n = size(A)
    indx = collect(1:m)
    sindx = zeros(eltype(A), m)
    nindx = zeros(eltype(A), m)
    col = zeros(eltype(A), m)
    for c in n:-1:1
        for i in 1:m
            col[i] = A[indx[i], c]
        end
        # Sorting a column vector is much faster than sorting a column matrix
        #sindx = sortperm(col, alg=QuickSort);
        sortperm!(sindx, col, alg = QuickSort) # available for 0.4, slightly faster
        for i in 1:m
            nindx[i] = indx[sindx[i]]
        end
        for i in 1:m
            indx[i] = nindx[i]
        end
    end
    return indx
end

function _countrepeats(A)
    m, n = size(A)
    occurrences = zeros(eltype(A), m)
    occurrences[1] = 1
    for i in 2:m
        if A[i, :] == A[i - 1, :]
            occurrences[i] = occurrences[i - 1] + 1
        else
            occurrences[i] = 1
        end
    end
    return occurrences
end

"""
    ir_skeleton(ir::IncRel, name = "skt")

Compute the "skeleton" of an incidence relation.

Compute a new incidence relation from an existing incidence relation to extract
the skeleton.

# Returns
Incidence relation for the skeleton of the input incidence relation. The left
shape collection consists of facets (shapes of manifold dimension one less than
the manifold dimension of the shapes themselves). The right shape collection is
the same as for the input.
"""
function ir_skeleton(ir::IncRel, name = "skt")
    @_check (manifdim(ir.right) == 0)
    @_check (manifdim(ir.left) > 0)
    inttype = eltype(ir._v[1])
    c = _asmatrix(ir, inttype) # incidence as a 2D array
    # construct a 2D array of the hyperface incidences
    hfc = c[:, facetconnectivity(ir.left, 1)]
    for i in 2:nfacets(shapedesc(ir.left))
        hfc = vcat(hfc, c[:, facetconnectivity(ir.left, i)])
    end
    # make a working copy
    s2hfc = deepcopy(hfc)
    _mysortdim2!(s2hfc) # sort columnwise
    idx = _mysortrowsperm(s2hfc) # find row permutation that sorts the rows
    # sort the rows of both the original hyper face array and the column-sorted array
    shfc = hfc[idx, :]
    s2hfc = s2hfc[idx, :]
    rep = _countrepeats(s2hfc) # find the repeats of the rows
    # identify boundary facets: the hyperface is boundary if it has no repeats
    isunq = [(rep[i] == 1) && (rep[i + 1] == 1) for i in 1:(length(rep) - 1)]
    push!(isunq, (rep[end] == 1)) # and the last one
    # unique rows are obtained by ignoring the repeats
    unq = findall(a -> a == 1, rep)
    unqhfc = shfc[unq, :] # unique hyper faces
    sir = IncRel(ShapeColl(facetdesc(ir.left), size(unqhfc, 1)), ir.right, unqhfc, name)
    sir.left.attributes["isboundary"] = VecAttrib(isunq[unq]) # store the is-boundary attribute
    return sir
end

"""
    ir_boundary(ir::IncRel, name = "bdr")

Compute the incidence relation for the boundary of the incidence relation on input.

The `skeleton` function.
"""
function ir_boundary(ir::IncRel, name = "bdr")
    sir = ir_skeleton(ir)
    isboundary = sir.left.attributes["isboundary"]
    ind = [i for i in 1:length(isboundary) if isboundary[i]]
    lft = ShapeColl(shapedesc(sir.left), length(ind), "facets")
    return IncRel(lft, sir.right, sir._v[ind], name)
end

function _selectrepeating(v, nrepeats)
    m = length(v)
    occurrences = zeros(eltype(v), m)
    sort!(v)
    occurrences[1] = 1
    for i in 2:m
        if v[i] == v[i - 1]
            occurrences[i] = occurrences[i - 1] + 1
        else
            occurrences[i] = 1
        end
    end
    repeats = [o == nrepeats for o in occurrences]
    return v[repeats]
end

"""
    ir_bbyfacets(ir::IncRel, fir::IncRel, tfir::IncRel, name = "bbf")

Compute the "bounded by facets" incidence relation `(d, d-1)` for `d`-dimensional shapes.

In other words, this is the incidence relation between shapes and the shapes
that bound these shapes (facets). For tetrahedra as the shapes, the incidence
relation lists the numbers of the faces that bound each individual tetrahedron.
The resulting left shape is of the same shape description as in the `ir`.

# Arguments
- `ir` = incidence relation `(d, 0)`,
- `fir` = incidence relation `(d-1, 0)`,
- `tfir` = transpose of the above, incidence relation `(0, d-1)`.

# Returns
Incidence relation for the bounded-by relationship between shape collections.
The left shape collection is the same as for the `ir`, the right shape
collection is for the facets (shapes of manifold dimension one less than the
manifold dimension of the shapes themselves).

!!! note
    The numbers of the facets are signed: positive when the facet bounds the shape
    in the sense in which it is defined by the shape as oriented with an outer
    normal; negative otherwise. The sense is defined by the numbering of the
    1st-order vertices of the facet shape.
"""
function ir_bbyfacets(ir::IncRel, fir::IncRel, tfir::IncRel, name = "bbf")
    @_check (manifdim(ir.right) == 0)
    @_check (manifdim(fir.right) == 0)
    @_check (manifdim(tfir.left) == 0)
    @_check manifdim(ir.left) == manifdim(tfir.right) + 1
    @_check manifdim(tfir.right) == manifdim(fir.left)
    inttype = eltype(ir._v[1])
    n1st = n1storderv(fir.left.shapedesc)
    nshif = nshifts(fir.left.shapedesc)
    # Sweep through the relations of (d, 0,) and use the (0, d-1) tfir
    _c = fill(inttype(0), nrelations(ir), nfacets(ir.left))
    for i in 1:nrelations(ir) # Sweep through the relations of (d, 0)
        sv = ir[i]
        c = inttype[] # this will be the list of facets at the vertices of this entity
        for j in 1:nentities(ir, i) # for all vertices
            fv = ir[i, j]
            for k in 1:nentities(tfir, fv)
                push!(c, tfir[fv, k])
            end # k
        end
        c = _selectrepeating(c, nvertices(tfir.right)) # keep the repeats
        @_check length(c) == nfacets(ir.left)
        # Now figure out the orientations
        for k in 1:nfacets(ir.left)
            fc = sv[facetconnectivity(ir.left, k)]
            sfc = sort(fc)
            for j in 1:length(c)
                oc = fir[c[j]]
                if sfc == sort(oc)
                    sgn = _sense(fc[1:n1st], oc, nshif)
                    _c[i, k] = sgn * c[j]
                    break
                end
            end # j
        end
    end
    bfacets = ShapeColl(fir.left.shapedesc, nshapes(tfir.right))
    return IncRel(ir.left, bfacets, _c, name)
end

"""
    ir_bbyfacets(ir::IncRel, fir::IncRel, name = "bbf")

Compute the "bounded by facets" incidence relation `(d, d-1)` for `d`-dimensional shapes.

Convenience function where the transpose of the incidence relation on the right
is computed on the fly.

# See also: [`ir_bbyfacets(ir::IncRel, fir::IncRel, tfir::IncRel)`](@ref)
"""
function ir_bbyfacets(ir::IncRel, fir::IncRel, name = "bbf")
    return ir_bbyfacets(ir, fir, ir_transpose(fir), name)
end

"""
    ir_bbyridges(ir::IncRel, eir::IncRel, teir::IncRel, name = "bbr")

Compute the "bounded by ridges" incidence relation `(d, d-2)` for `d`-dimensional shapes.

In other words, this is the incidence between shapes and the shapes that "bound"
the boundaries of these shapes (i. e. ridges). For tetrahedra as the shapes, the
incidence relation lists the numbers of the edges that "bound" each individual
tetrahedron. The resulting shape is of the same shape description as the
`shapes` on input.

# Returns
The incidence relation for the "bounded by" relationship between shape
collections. The left shape collection is the same as for the `ir`, the right
shape collection is for the ridges (shapes of manifold dimension two less than
the manifold dimension of the shapes themselves).

!!! note
    The numbers of the ridges are signed: positive when the ridge bounds the shape
    in the sense in which it is defined by the shape as oriented with an outer
    normal; negative otherwise. The sense is defined by the numbering of the
    1st-order vertices of the ridge shape.
"""
function ir_bbyridges(ir::IncRel, eir::IncRel, teir::IncRel, name = "bbr")
    @_check (manifdim(ir.right) == 0) && (manifdim(teir.left) == 0)
    @_check manifdim(ir.left) == manifdim(teir.right) + 2
    @_check manifdim(teir.right) == manifdim(eir.left)
    inttype = eltype(ir._v[1])
    n1st = n1storderv(eir.left.shapedesc)
    nshif = nshifts(eir.left.shapedesc)
    # Sweep through the relations of (d, 0), and use the (0, d-2) teir
    _c = fill(inttype(0), nrelations(ir), nridges(ir.left))
    for i in 1:nrelations(ir) # Sweep through the relations of (d, 0)
        sv = ir[i]
        c = inttype[] # this will be the list of facets at the vertices of this entity
        for j in 1:nentities(ir, i) # for all vertices
            fv = ir[i, j]
            for k in 1:nentities(teir, fv)
                push!(c, teir[fv, k])
            end # k
        end
        c = _selectrepeating(c, nvertices(teir.right)) # keep the repeats
        @_check length(c) == nridges(ir.left)
        # Now figure out the orientations
        for k in 1:nridges(ir.left)
            fc = sv[ridgeconnectivity(ir.left, k)]
            sfc = sort(fc)
            for j in 1:length(c)
                oc = eir[c[j]]
                if sfc == sort(oc)
                    sgn = _sense(fc[1:n1st], oc, nshif)
                    _c[i, k] = sgn * c[j]
                    break
                end
            end # j
        end
    end
    bridges = ShapeColl(eir.left.shapedesc, nshapes(teir.right))
    return IncRel(ir.left, bridges, _c, name)
end

"""
    ir_bbyridges(ir::IncRel, eir::IncRel, name = "bbr")

Compute the "bounded-by-ridges" incidence relation `(d, d-2)`.

Convenience function where the transpose of the incidence relation on the right
is computed on the fly.

# See also: [`ir_bbyridges(ir::IncRel, eir::IncRel, efir::IncRel)`](@ref)
"""
function ir_bbyridges(ir::IncRel, eir::IncRel, name = "bbr")
    return ir_bbyridges(ir, eir, ir_transpose(eir), name)
end

"""
    ir_identity(ir::IncRel, side = :left)

Compute the "identity" incidence relation `(d, d)`.

Either the left (side = :left) or the right (side = :right) shape is mapped
to itself.
"""
function ir_identity(ir::IncRel, side = :left)
    if side == :left
        return IncRel(
            ir.left,
            ir.left,
            [SVector{1,Int64}([idx]) for idx in 1:nshapes(ir.left)],
        )
    else
        return IncRel(
            ir.right,
            ir.right,
            [SVector{1,Int64}([idx]) for idx in 1:nshapes(ir.right)],
        )
    end
end

"""
    ir_subset(ir::IncRel, list)

Compute the incidence relation representing a subset of the shape collection
on the left.
"""
function ir_subset(ir::IncRel, list)
    ssl = ShapeColl(ir.left.shapedesc, length(list))
    return IncRel(ssl, ir.right, ir._v[list], ir.left.name * "subset")
end

"""
    Base.summary(ir::IncRel)

Form a brief summary of the incidence relation.
"""
function Base.summary(ir::IncRel)
    return "$(ir.name): " * summary(ir.left) * ", " * summary(ir.right)
end


text:
import FEM.stiffness
import FEM.Bmatrix
import FEM.creategps
import FEM.createinterp
import FEM.createstorage
import FEM.create_matstat
import FEM.Vertex3
import FEM.Vertex4
facts("FEM.Element") do

FEM.lintrigmod()
using FEM.LinTrigMod

FEM.linquadmod()
using FEM.LinQuadMod

FEM.linearisotropicmod()
using FEM.LinearIsotropicMod

context("FEM.Element.LinTrig") do

    nodes = [FENode2(1, [0, 0]), FENode2(2, [1, 1]), FENode2(3, [1, 2])]

    # Create the gps, interpolator, storage...

    gps = creategps(LinTrig)
    interp = createinterp(LinTrig)
    storage = createstorage(LinTrig)
    mat = LinearIsotropic(200e9, 0.3)
    ms = create_matstat(typeof(mat))
    elem = LinTrig(Vertex3(1, 2, 3), 1, interp, storage, gps, ms)

    Ke = stiffness(elem, nodes, mat, :dof_vals)

    # Calculated with plante([0.0, 1.0, 1.0], [0.0, 1.0, 2.0], [2, 1], hooke(2, 200e9, 0.3))
    Ke_calfem =  1.0e+11 * [
    [ 1.346153846153846                   0  -2.692307692307692   0.576923076923077   1.346153846153846  -0.576923076923077];
    [                 0   0.384615384615385   0.384615384615385  -0.769230769230769  -0.384615384615385   0.384615384615385];
    [-2.692307692307692   0.384615384615385   5.769230769230769  -1.923076923076923  -3.076923076923077   1.538461538461538];
    [ 0.576923076923077  -0.769230769230769  -1.923076923076923   2.884615384615385   1.346153846153846  -2.115384615384615];
    [ 1.346153846153846  -0.384615384615385  -3.076923076923077   1.346153846153846   1.730769230769231  -0.961538461538461];
    [-0.576923076923077   0.384615384615385   1.538461538461538  -2.115384615384615  -0.961538461538461   1.730769230769231]]

    @fact norm(Ke - Ke_calfem) / norm(Ke) => roughly(0.0)

end # context

context("FEM.Element.LinQuad") do

    nodes = [FENode2(1, [0.0, 0.0]), FENode2(2, [1.0, 0.0]), FENode2(3, [1.0, 2.5]), FENode2(4, [0.0, 1.5])]

    gps = creategps(LinQuad)
    interp = createinterp(LinQuad)
    storage = createstorage(LinQuad)
    mat = LinearIsotropic(200e9, 0.3)
    ms = create_matstat(typeof(mat))
    elem = LinQuad(Vertex4(1, 2, 3, 4), 1, interp, storage, gps, ms)

    Ke = stiffness(elem, nodes, mat, :dof_vals)

    # Calculated with plani4e([0.0, 1.0, 1.0, 0.0], [0.0, 0.0, 2.5, 1.5], [2, 1, 2], hooke(2, 200e9, 0.3))
    Ke_calfem = 1.0e+11 * [
    [ 1.665302782324059   0.296644844517185  -1.653027823240589   0.014320785597381  -0.366202945990180  -0.398936170212766   0.353927986906710   0.087970540098200];
    [ 0.296644844517185   0.949263502454992  -0.177986906710311  -0.261865793780688  -0.398936170212766  -0.315057283142389   0.280278232405892  -0.372340425531915];
    [-1.653027823240589  -0.177986906710311   2.491816693944353  -0.623977086743044   0.873567921440262   0.239361702127660  -1.712356792144026   0.562602291325696];
    [ 0.014320785597381  -0.261865793780687  -0.623977086743044   1.080196399345336   0.047054009819967  -0.118657937806874   0.562602291325696  -0.699672667757774];
    [-0.366202945990180  -0.398936170212766   0.873567921440262   0.047054009819967   1.145662847790507   0.337561374795417  -1.653027823240589   0.014320785597381];
    [-0.398936170212766  -0.315057283142389   0.239361702127660  -0.118657937806874   0.337561374795417   0.695581014729951  -0.177986906710311  -0.261865793780688];
    [ 0.353927986906710   0.280278232405892  -1.712356792144026   0.562602291325696  -1.653027823240589  -0.177986906710311   3.011456628477905  -0.664893617021277];
    [ 0.087970540098200  -0.372340425531915   0.562602291325696  -0.699672667757774   0.014320785597381  -0.261865793780687  -0.664893617021277   1.333878887070377]]

    @fact norm(Ke - Ke_calfem) / norm(Ke) => roughly(0.0, 10.0^(-15))

end # context

end # facts

ext: jl


text:
using Extremes, DataFrames, Gadfly, Test, Distributions

data = load("portpirie")

x = Variable("x", randn(100))
μ = 10 .+ x.value
σ = 1.0
ξ = .1
pd = GeneralizedExtremeValue.(μ, σ, ξ)
y = rand.(pd)

fm = gevfit(data, :SeaLevel)
fmns = gevfit(y, locationcov=[x])

probplot(fm)
probplot(fmns)

qqplot(fm)
qqplot(fmns)

returnlevelplot(fm)
returnlevelplot(fmns)

histplot(fm)
histplot(fmns)

diagnosticplots(fm)
diagnosticplots(fmns)

fm = gevfitbayes(data, :SeaLevel)
fmns = gevfitbayes(y, locationcov=[x])

probplot(fm)
probplot(fmns)

qqplot(fm)
qqplot(fmns)

returnlevelplot(fm)
returnlevelplot(fmns)

histplot(fm)
histplot(fmns)

diagnosticplots(fm)
diagnosticplots(fmns)

ys = rand(GeneralizedPareto(1,.1),100)

x = Variable("x", randn(100)/3)
ϕ = x.value
σ = exp.(ϕ)
ξ = .1
pd = GeneralizedPareto.(σ, ξ)
y = rand.(pd)

fm = gpfit(ys)
fmns = gpfit(y, logscalecov=[x])

probplot(fm)
probplot(fmns)

qqplot(fm)
qqplot(fmns)

returnlevelplot(fm)
returnlevelplot(fmns)

histplot(fm)
histplot(fmns)

diagnosticplots(fm)
diagnosticplots(fmns)

fm = gpfitbayes(ys)
fmns = gpfitbayes(y, logscalecov=[x])

probplot(fm)
probplot(fmns)

qqplot(fm)
qqplot(fmns)

returnlevelplot(fm)
returnlevelplot(fmns)

histplot(fm)
histplot(fmns)

diagnosticplots(fm)
diagnosticplots(fmns)


text:
###### variogram estimator ######

"""
Multifractal estimator.

# Args
- estimator: function of estimator having the interface `estimator(X, p, args...; kwargs...)`. Here `X` is the input array, `p` is the power of moment.
"""
function multifractal_estim(estimator::Function, X::AbstractArray, pows::AbstractVector{<:Real}, args...; kwargs...)
    @assert all(pows .> 0) "Powers of moment must be positive"

    res = [estimator(X, p, args...; kwargs...) for p in pows]
    ws = StatsBase.weights([x.residual for x in res].^ -1)
    Hs = [x.hurst for x in res]
    σs = [x.σ for x in res]
    Rs = [x.residual for x in res]

    # return (hurst=mean(Hs, ws), σ=mean(σs, ws), residual=mean(Rs, ws))
    # return (hurst=median(Hs, ws), σ=median(σs, ws), residual=median(Rs, ws))
    return (hurst=median(Hs), σ=median(σs), residual=median(Rs))
end

"""
Variance function of the empirical variogram
"""
function variogram_variance(H::Real, lags::AbstractVector{<:Integer}, N::Integer)
    V = zeros(Float64, length(lags))

    for (j, d) in enumerate(lags)
        proc = FractionalGaussianNoise(H, d)
        V[j] = autocov(proc, 0)^2 + sum((1-i/N) * autocov(proc, i)^2 for i=0:N-1)
    end
    return 4/N * V
end

"""
    powlaw_estim(X, lags; p=2., method=:optim)

Power-law estimator for Hurst exponent and volatility.

# Args
- X: matrix of fGn, each row is a fGn of some time lag and each column is an observation.
- lags: time lags (increment step) used to compute each component of `X`
- p: power of the moment
- methods: method of estimation: {:optim, :lm}

# Returns
- hurst, σ: estimation of Hurst and volatility, as well as an object of optimizer
- (xp, yp): vectors of regression
- opm: optimizer

# Notes
- `X` is computed from fBm by taking finite differences. The second dimension corresponds to time. Example, let `W` be a fBm sample path then the following command computes `X`:
```julia
julia> lags = 2:10
julia> X = transpose(lagdiff(W, lags, mode=:causal))
```
- `p=1` is robust against quantization error.
"""
function powlaw_estim(X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; pow::Real=2., method::Symbol=:lm, reweight::Bool=true)
    # remove columns containing NaN
    idx = findall(vec(.!any(isnan.(X), dims=1)))
    X = X[:,idx] # view(X,:,idx)

    if length(X) == 0
        return (hurst=NaN, σ=NaN, η=NaN, residual=NaN, vars=nothing, optimizer=nothing)
    else
        @assert length(lags) == size(X,1) > 1  "Dimension mismatch."
        @assert all(lags .>= 1)  "Lags must be larger than or equal to 1."
        @assert pow>0  "Moment must be positive."

        # explanatory and observation vectors
        xp = pow * log.(lags)
        μX = mean(X, dims=2)
        yp = vec(log.(mean((abs.(X.-μX)).^pow, dims=2)))
        # yp = vec(log.(mean((abs.(X)).^pow, dims=2)))
        # hurst, η, res, err = NaN, NaN, NaN, (NaN, NaN)

        # old version with weights
        # # polynomial order of the weight for samples, if 0 the uniform weight is used
        # kt::Integer = 0  # non-zero value puts more weight on most recent samples (i.e. those at large column numbers).
        # wt = StatsBase.weights(causal_weight(size(X,2), kt))
        # μX = mean(X, wt, 2)
        # yp = vec(log.(mean((abs.(X.-μX)).^p, wt, 2)))
        # # yp = vec(log.(mean((abs.(X)).^p, wt, 2)))  # <- this gives lower SNR
        # xp = p * log.(lags)
        # # weight for scales
        # ks::Integer = 0  # hard-coded: polynomial order of weight for scales, if 0 the uniform weight is used
        # ws = StatsBase.weights(poly_weight(length(yp), ks))
        # yc = yp .- mean(yp, ws)
        # xc = xp .- mean(xp, ws)
        # func = h -> 1/2 * sum(ws .* (yc - h*xc).^2)
        # # func = h -> 1/2 * sum(ws .* abs.(yc - h*xc))

        # compute the weighting vector:
        # Run first an estimation of Hurst by linear regression, then use this estimate to compute the weighting vector.
        # Note that the reweighting scheme is exact only for pow=2.
        dg = DataFrames.DataFrame(xvar=xp, yvar=yp)
        opm = GLM.lm(@GLM.formula(yvar~xvar), dg)
        η, hurst = GLM.coef(opm)  # intercept and slope
        res = GLM.deviance(opm)  # residual
        # err = try  # std error of estimates
        #     GLM.stderror(opm)
        # catch
        #     (NaN, NaN)
        # end

        ws =  if reweight
            variogram_variance(0 < hurst < 1 ? hurst : 0.5, lags, size(X,2)) .^ (-1)
            # variogram_variance(0.5, lags, size(X,2)) .^ -1  # in practice this is good enough?
        else
            ones(length(lags))  # uniform weight
        end
        ws ./= sum(ws)

        # estimation of H and η
        if method == :optim
            yc = yp .- mean(yp)
            xc = xp .- mean(xp)
            func = h -> sum(ws .* (yc - h*xc).^2)
            # func = h -> sum(ws .* abs.(yc - h*xc).^qnorm)
            # Gradient-free constrained optimization
            # ɛ = 1e-2  # search hurst in the interval [ɛ, 1-ɛ]
            opm = Optim.optimize(func, 0., 1., Optim.Brent())
            # # Gradient-based optimization
            # optimizer = Optim.GradientDescent()  # e.g. Optim.BFGS(), Optim.GradientDescent()
            # opm = Optim.optimize(func, ε, 1-ε, [0.5], Optim.Fminbox(optimizer))
            hurst = Optim.minimizer(opm)[1]
            η = mean(yp - hurst*xp)
            # relative residual
            res = sqrt(opm.minimum / length(xp) / var(yp))
        elseif method == :lm  # using GLM package
            dg = DataFrames.DataFrame(xvar=xp, yvar=yp)
            opm = GLM.glm(@GLM.formula(yvar~xvar), dg, GLM.Normal(), GLM.IdentityLink(), wts=sqrt.(ws))
            # opm = GLM.lm(@GLM.formula(yvar~xvar), dg)
            η, hurst = GLM.coef(opm)
            # GLM.deviance is by definition the RSS
            res = sqrt(GLM.deviance(opm) / var(yp)) # / length(xp)

            # # or equivalently, by manual inversion
            # Ap = hcat(xp, ones(length(xp))) # design matrix
            # hurst, η = Ap \ yp
        else
            error("Unknown method $(method).")
        end

        cp = normal_moment_factor(pow)  # constant factor depending on p
        σ = exp((η-log(cp))/pow)

        # return hurst, σ, (xp, yp)
        return (hurst=hurst, σ=σ, η=η, residual=res, vars=(xp,yp), optimizer=opm)
    end
end

"""
    powlaw_estim(X::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; kwargs...)

# Args
- X: sample path of fBm, e.g. the log-price
- lags: time lags used to compute finite differences
"""
function powlaw_estim(X::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; mode::Symbol=:causal, kwargs...)
    dX = transpose(lagdiff(X, lags; mode=mode))  # take transpose s.t. each column is an observation
    return powlaw_estim(dX, lags; kwargs...)
end

const variogram_estim = powlaw_estim

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
using Test
using DataTypesBasic

@test convert(Const, []) == Const(nothing)
@test_throws AssertionError convert(Const, [1,2,3])

@test convert(Vector, Identity("hello")) == ["hello"]
@test convert(Vector, Const(43)) == []
@test convert(Vector, @ContextManager cont -> cont(2)) == [2]
@test_throws MethodError convert(Vector, nothing)

@test convert(Identity, @ContextManager cont -> cont(2)) == Identity(2)
@test_throws MethodError convert(Identity, [1,2,3])

@test run(convert(ContextManager, Identity(4))) == 4
@test_throws MethodError convert(ContextManager, nothing)
@test_throws MethodError convert(ContextManager, Const(:error))
@test_throws MethodError convert(ContextManager, [1,2,3])

@test convert(Either{String, Float32}, Identity(1)) == Identity(Float32(1))
@test convert(Either{String, Int}, Identity(1)) == Identity(1)
@test_throws MethodError convert(Either{String, Int}, Identity("hi"))

@test convert(Either{Float32, Int}, Const(1)) == Const(Float32(1))
@test convert(Either{Int, Int}, Const(1)) == Const(Float32(1))
@test_throws MethodError convert(Either{String, Int}, Const(1))

@test convert(Either{String}, Identity(1)) == Identity(1)
@test convert(Either{String}, Const("hi")) == Const("hi")
@test convert(Either{Float32}, Const(1)) == Const(Float32(1))
@test_throws MethodError convert(Either{Float32}, Const("hi"))

@test convert(Either{<:Any, String}, Const(1)) == Const(1)
@test convert(Either{<:Any, String}, Identity("hi")) == Identity("hi")
@test convert(Either{<:Any, Float32}, Identity(1)) == Identity(Float32(1))
@test_throws MethodError convert(Either{<:Any, Float32}, Identity("hi"))


text:
module InsertInBlocksTests

using Test
using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.Fields
using Gridap.TensorValues
using Gridap.Fields: MockField, MockBasis

l = 10
np = 3
ndofs1 = 4
ndofs2 = 6

p = Point(1,2)
x = fill(p,np)
z = 2.0

v = VectorValue(3.0,1.5)
w = VectorValue(3.4,3.5)
a = MockBasis{2}(v,ndofs2)
b = MockBasis{2}(w,ndofs1)

al = Fill(a,l)
atl = trialize_array_of_bases(al)
bl = Fill(b,l)
xl = Fill(x,l)

ax1 = Fill((Base.OneTo(ndofs1),),l)
ax2 = Fill((Base.OneTo(ndofs2),),l)
aBl = insert_array_of_bases_in_block(2,al,ax1,ax2)
aBl_x = evaluate(aBl,xl)
@test isa(aBl,VectorOfBlockBasisCoo)
@test isa(aBl_x,VectorOfBlockArrayCoo)

ax1 = Fill((Base.OneTo(ndofs1),),l)
ax2 = Fill((Base.OneTo(ndofs2),),l)
aBl = insert_array_of_bases_in_block(3,al,ax1,ax2,ax2)
aBl_x = evaluate(aBl,xl)
@test isa(aBl,VectorOfBlockBasisCoo)
@test isa(aBl_x,VectorOfBlockArrayCoo)

bBl = insert_array_of_bases_in_block(1,bl,ax1,ax2)
bBl_x = evaluate(bBl,xl)
@test isa(bBl,VectorOfBlockBasisCoo)
@test isa(bBl_x,VectorOfBlockArrayCoo)

ax1 = Fill((Base.OneTo(1),Base.OneTo(ndofs1)),l)
ax2 = Fill((Base.OneTo(1),Base.OneTo(ndofs2)),l)
atBl = insert_array_of_bases_in_block(2,atl,ax1,ax2)
atBl_x = evaluate(atBl,xl)
@test isa(atBl,VectorOfBlockBasisCoo)
@test isa(atBl_x,VectorOfBlockArrayCoo)

ax1 = Fill((Base.OneTo(1),Base.OneTo(ndofs1)),l)
ax2 = Fill((Base.OneTo(1),Base.OneTo(ndofs2)),l)
atBl = insert_array_of_bases_in_block(3,atl,ax1,ax2,ax2)
atBl_x = evaluate(atBl,xl)
@test isa(atBl,VectorOfBlockBasisCoo)
@test isa(atBl_x,VectorOfBlockArrayCoo)

ax1 = aBl.axes
ax2 = aBl.axes
aSl = insert_array_of_bases_in_block(2,aBl,ax1,ax2)
@test isa(aSl,VectorOfBlockBasisCoo)
aSl_x = evaluate(aSl,xl)
@test isa(aSl_x,VectorOfBlockArrayCoo)

ax1 = atBl.axes
ax2 = atBl.axes
atSl = insert_array_of_bases_in_block(1,atBl,ax1,ax2)
@test isa(atSl,VectorOfBlockBasisCoo)
atSl_x = evaluate(atSl,xl)
@test isa(atSl_x,VectorOfBlockArrayCoo)

blocks = (al,)
blockids = [(2,)]
ranges_i = (blockedrange([ndofs1,ndofs2]),)
ranges = Fill(ranges_i,l)
aBl = VectorOfBlockBasisCoo(blocks,blockids,ranges)

blocks = (bl,)
blockids = [(1,)]
bBl = VectorOfBlockBasisCoo(blocks,blockids,ranges)

blocks = (atl,)
blockids = [(1,2)]
ranges_i = (blockedrange([1]),blockedrange([ndofs1,ndofs2]))
ranges = Fill(ranges_i,l)
atBl = VectorOfBlockBasisCoo(blocks,blockids,ranges)

aBl_x = evaluate(aBl,xl)
@test isa(aBl_x,VectorOfBlockArrayCoo)
@test isa(aBl_x.axes,Fill)
@test aBl_x.axes[1] == (blockedrange([np]),blockedrange([ndofs1,ndofs2]))

bBl_x = evaluate(bBl,xl)
@test isa(bBl_x,VectorOfBlockArrayCoo)

atBl_x = evaluate(atBl,xl)
@test isa(atBl_x,VectorOfBlockArrayCoo)
@test isa(atBl_x.axes,Fill)
@test atBl_x.axes[1] == (blockedrange([np]),blockedrange([1]),blockedrange([ndofs1,ndofs2]))

atBl = trialize_array_of_bases(aBl)
atBl_x = evaluate(atBl,xl)
@test isa(atBl_x,VectorOfBlockArrayCoo)

btBl = trialize_array_of_bases(bBl)
btBl_x = evaluate(btBl,xl)
@test isa(btBl_x,VectorOfBlockArrayCoo)

∇aBl = ∇(aBl)
@test isa(∇aBl,VectorOfBlockBasisCoo)
∇aBl_x = evaluate(∇aBl,xl)
@test isa(∇aBl_x,VectorOfBlockArrayCoo)

cl = operate_arrays_of_fields(-,aBl,bBl)
cl_x = evaluate(cl,xl)
@test isa(cl_x,VectorOfBlockArrayCoo)
∇cl = ∇(cl)
∇cl_x = evaluate(∇cl,xl)
@test isa(cl_x,VectorOfBlockArrayCoo)

cl = operate_arrays_of_fields(⋅,atBl,cl)
cl_x = evaluate(cl,xl)
@test isa(cl_x,VectorOfBlockArrayCoo)
@test is_zero_block(cl_x,Block(1,1,1))
@test is_nonzero_block(cl_x,Block(1,1,2))
@test is_zero_block(cl_x,Block(1,2,1))
@test is_nonzero_block(cl_x,Block(1,2,2))

ids = [2,3,5,1,2]
cj = reindex(aBl,ids)
xj = reindex(xl,ids)
@test isa(cj,VectorOfBlockBasisCoo)
cj_x = evaluate(cj,xj)
@test isa(cj_x.axes,Fill)
@test cj_x.axes[1] == (blockedrange([np]),blockedrange([ndofs1,ndofs2]))
@test length(cj_x) == length(ids)
@test length(cj) == length(ids)

v = VectorValue(1.0,2.0,3.0)
f = MockField{2}(v)
fl = Fill(f,l)

cl = compose_field_arrays(aBl,fl)
@test isa(cl,VectorOfBlockBasisCoo)

#cache = array_cache(cl_x)
#using BenchmarkTools
#@btime getindex!($cache,$cl_x,3)

end # module


text:
# collection of calls

# /**
#  * rsvg_handle_new_from_stream_sync:
#  * @input_stream: a #GInputStream
#  * @base_file: (allow-none): a #GFile, or %NULL
#  * @flags: flags from #RsvgHandleFlags
#  * @cancellable: (allow-none): a #GCancellable, or %NULL
#  * @error: (allow-none): a location to store a #GError, or %NULL
#  *
#  * Creates a new #RsvgHandle for @stream.
#  *
#  * If @cancellable is not %NULL, then the operation can be cancelled by
#  * triggering the cancellable object from another thread. If the
#  * operation was cancelled, the error %G_IO_ERROR_CANCELLED will be
#  * returned.
#  *
#  * Returns: a new #RsvgHandle on success, or %NULL with @error filled in
#  *
#  * Since: 2.32
#  */
# RsvgHandle *
# rsvg_handle_new_from_stream_sync (GInputStream   *input_stream,
#                                   GFile          *base_file,
#                                   RsvgHandleFlags flags,
#                                   GCancellable    *cancellable,
#                                   GError         **error)
# {

"""
handle_get_dimensions(handle::RsvgHandle, dimension_data::RsvgDimensionData)
"""
function handle_get_dimensions(handle::RsvgHandle, dimension_data::RsvgDimensionData)
    ccall((:rsvg_handle_get_dimensions,librsvg), Nothing,
        (RsvgHandle,Ptr{RsvgDimensionData}), handle, Ref(dimension_data))
end

"""
dimension_data = handle_get_dimensions(handle::RsvgHandle)
"""
function handle_get_dimensions(handle::RsvgHandle)
    d = Rsvg.RsvgDimensionData(1,1,1,1);
    handle_get_dimensions(handle, d);
    d
end

"""
set_default_dpi(dpi::Float64)
"""
function set_default_dpi(dpi::Float64)
    ccall((:rsvg_set_default_dpi, librsvg), Nothing,
            (Float64,), dpi)
end
"""
handle_set_dpi(handle::RsvgHandle, dpi::Float64)
"""
function handle_set_dpi(handle::RsvgHandle, dpi::Float64)
    ccall((:rsvg_handle_set_dpi, librsvg), Nothing,
            (RsvgHandle,Float64), handle, dpi)
end

"""
handle_render_cairo(cr::CairoContext, handle::RsvgHandle)
"""
function handle_render_cairo(cr::CairoContext, handle::RsvgHandle)
	ccall((:rsvg_handle_render_cairo, librsvg), Bool,
        (RsvgHandle,Ptr{Nothing}), handle, cr.ptr)
end

"""
handle_new_from_file(filename::AbstractString,error::GError)
"""
function handle_new_from_file(filename::AbstractString,error::GError)
    ptr = ccall((:rsvg_handle_new_from_file, librsvg), Ptr{Nothing},
                (Cstring,GError), string(filename), error)
    RsvgHandle(ptr)
end

"""
handle_new_from_file(filename::AbstractString)
"""
handle_new_from_file(filename::AbstractString) = handle_new_from_file(filename::AbstractString,GError(0,0,0))

"""
handle_new_from_data(data::AbstractString,error::GError)
"""
function handle_new_from_data(data::AbstractString,error::GError)
    ptr = ccall((:rsvg_handle_new_from_data, librsvg), Ptr{Nothing},
                (Ptr{UInt8},UInt32,GError), string(data), sizeof(data),error)
    RsvgHandle(ptr)
end

"""
handle_new_from_data(data::AbstractString)
"""
function handle_new_from_data(data::AbstractString)

    perr = Ref{Ptr{GError}}(C_NULL)

    g_input_stream = Rsvg.glib_memory_input_stream_new_from_data(data)

    ptr = ccall((:rsvg_handle_new_from_stream_sync, librsvg), Ptr{Nothing},
        (GInputStream,Ptr{Nothing},UInt32,Ptr{Nothing},Ptr{Ptr{GError}}),
        g_input_stream,C_NULL,0,C_NULL,perr)

    if ptr == C_NULL
        err = unsafe_load(perr[])
        message = unsafe_string(err.message)
        ccall((:g_error_free, librsvg), Cvoid, (Ptr{GError},), perr[])
        error(message)
    end

    RsvgHandle(ptr)
end

"""
handle_free(handle::RsvgHandle)
"""
function handle_free(handle::RsvgHandle)
    ccall((:rsvg_handle_free,librsvg), Nothing, (RsvgHandle,), handle)
end

"""
handle_free(handle::GInputStream)
"""
function handle_free(handle::GInputStream)
    ccall((:g_input_stream_close,libgio), Nothing,
        (GInputStream,Ptr{Nothing},Ptr{Nothing}), handle,C_NULL,C_NULL)
end


text:
using ExpressWorkflowMaker
using Test

@testset "ExpressWorkflowMaker.jl" begin
    # Write your tests here.
end


text:
@userplot espplott

@recipe function f(h::espplott;
        n=15, stg=false, fq=false, VisCor=false,μ=1)
     Esp, x = h.args

   CapRange, Head, EspPower, η, Cone, Stages, Freq, CapCorrected, HeadCorrected, ηCorrected, PowerCorrected = EspPerformance(Esp,ns=n,stgs=stg, fqs=fq, VisCors=VisCor, μs=μ)

    xformatter := :plain
    legend := :topright
    layout := (3,1)
    link := :x
    size := (800,600)
    legendfontsize := 6
    xlim:=(0,Esp.AOF)
    if VisCor==true
        @series begin
            seriestype := :path
            linestyle := :dash
            linewidth := 1
            linecolor := :darkblue
            marker := :hexagon
            markercolor := :blue
            subplot :=1
            label := "Head Corrected"
            CapCorrected, HeadCorrected
        end

        @series begin
            seriestype := :path
            linestyle := :dash
            linewidth := 1
            linecolor := :darkblue
            marker := :hexagon
            markercolor := :black
            subplot :=2
            label := "Eff Corrected"
            CapCorrected, ηCorrected
        end
        @series begin
            seriestype := :path
            linestyle := :dash
            linewidth := 1
            linecolor := :darkgreen
            marker := :hexagon
            markercolor := :green
            subplot :=3
            label := "Power Corrected"
            CapCorrected, PowerCorrected
        end

    end
    @series begin
        seriestype := :path
        linewidth := 2
        linecolor := :blue
        subplot :=1
        label := "Head"
        CapRange, Head
    end

    @series begin
        seriestype := :scatter
        subplot :=1
        label :="Min/Max"
        Cone.caps, Cone.heads
    end

         @series begin
        seriestype := :vline
        linewidth := 1
        linestyle := :dash
        linecolor := :red
        subplot :=1
        label := "Efficiency Range"
        ylabel := "Head [ft]"
        title := " ESP Performance Curve \n $(Esp.Manufacter) Series $(Esp.Series) Model $(Esp.Model) \n at $(Freq) Hz and $(Stages) Stages"
        titlefontsize := 8
        Cone.caps
    end

    @series begin
        seriestype := :path
        linewidth := 1
        linecolor := :black
        subplot :=2
        label := "Efficiency"
         ylabel := "Efficiency [%]"
        CapRange, η
    end

 #
      @series begin
        seriestype := :path
        link := :x
        linewidth := 2
        linecolor := :darkgreen
        label:= "Horse Power"
        subplot :=3
        xlabel := "Liquid Rate [bbl/d]"
         ylabel := "Horse Power [hp]"
        CapRange,EspPower
    end

end

@userplot espmatch

@recipe function f(h::espmatch;
                            pip=false, n=15, rfq=30:5:60, stg=false, fq=false,
                            pbs=0,
                            ϵs=0.0006, Tss=80, ∇Ts=1, GORs=false, dis=2.99, Dns=50, tols=0.05, sgs=0.6)

        Esp, Pr, J, OilPVT, GasPVT, WaterPVT, Thp, PumpDepth, PerfDepth, Ql, Bsw=h.args

    CapRange,Head,CapMin,CapMax,HeadMin,HeadMax, CapCustom, HeadCustom=EspTornado(Esp; ns=n, rfqs=rfq, stgs=stg, fqs=fq)

        ΔD=PerfDepth-PumpDepth

        ρo=LinearInterpolation(OilPVT.P.*1, OilPVT.Rho.*1)
        ρw=LinearInterpolation(WaterPVT.P.*1, WaterPVT.Rho.*1)
        ρl= (ρo(Pr)*(1-Bsw)+ρw(Pr)*Bsw)
        ∇P=ρl*(0.433/62.4)
        if pip == false
          Pipr, Qipr, leg=OilInflow(Pr,J,pb=pbs,n=20)
          Pinter=LinearInterpolation(reverse(vec(Qipr)),reverse(vec(Pipr)))      # Interpolation accepts increasing order and vectors only
          pwf=Pinter(Ql)
          pip=pwf-∇P*ΔD
        end

        println("Pump Intake Pressure $(round(pip, digits=1)) psi")
        DepthRange, Qrange, Pd, legO, P=OilOutflowSen(PumpDepth,Thp,Bsw,OilPVT,GasPVT,WaterPVT, CapCustom;ϵ=ϵs, Ts=Tss, ∇T=∇Ts, GOR=GORs, di=dis,  Dn=Dns, tol=tols, sg=sgs)

        HeadMatch=(Pd'.-pip)./∇P

       xlabel := "Liquid Rate [bbl/d]"
       ylabel := "Head [ft]"
       xformatter := :plain

        @series begin
            seriestype := :path
            linecolor := :black
            linewidth := 1
            label := map(x->"$x Hz",rfq)
            CapRange,Head
        end

        @series begin
            seriestype := :path
            linecolor := :darkgreen
            linestyle := :dash
            linewidth := 3
            label := "At $fq Hz"
            CapCustom, HeadCustom
        end

        @series begin
            seriestype := :path
            linecolor := :red
            linewidth := 2
            linestyle := :dash
            label := "Min/Max"
            [CapMin, CapMax],[HeadMin, HeadMax]
        end

        @series begin
                seriestype := :path
                linecolor := :green
                linewidth := 2
                label := "ESP Match"
                CapCustom, HeadMatch
        end

    end


text:
using Plots, LaTeXStrings, Statistics

function Deposition_2DG(; len, tot_time, dep_rate, color_step)
    surf = [[] for i=1:len]
    for i in surf
        push!(i, [0, 0])
    end
    for n in 0:tot_time
        randsurf = rand(1:len, dep_rate)
        for i in randsurf
            maxlen = max(surf[sides(i, len)[1]][end][1], surf[sides(i, len)[2]][end][1], surf[i][end][1]+1)
            push!(surf[i], [maxlen, n%color_step])
        end
    end
    return Leveling(surf, len)
end

function sides(n, L)
    if n == L
        return n-1 , 1
    elseif n == 1
        return L , n+1
    else
        return n-1, n+1
    end
end

function Leveling(surf, len)
    lasth = []
    for i in surf
        append!(lasth, i[end][1])
    end

    surface = zeros(Int,(len, max(lasth...)+1))

    for i in 1:len
        for j in surf[i]
            surface[i, j[1]+1] = j[2]
        end
    end
    return transpose(hcat(surface))
end

Parameters = Dict(:len => 200,
                    :tot_time => 100,
                        :dep_rate => 1000,
                            :color_step => 10)

surf = Deposition_2DG(; Parameters...)

theme(:dark)
gr()

heatmap(surf, c = :solar, legend = false)
savefig("C:\\Users\\Yaghoub\\Documents\\GitHub\\Ballistic-Deposition\\Deposition\\Ballistic-Deposition-Kardar-Parisi-Zhang\\Fig\\deposition.png")

max_stars_repo_name: shahmari/Ballistic-Deposition

max_issues_repo_name: shahmari/Ballistic-Deposition

max_forks_repo_name: shahmari/Ballistic-Deposition


text:
#!/usr/local/bin/julia
# last tested Julia version: 1.6.1
# Complex Stuart-Landau Equation model on a 2D lattice
# FvW 03/2018

using Interpolations
using NPZ
using PyCall
using PyPlot
using Statistics
using VideoIO
@pyimport matplotlib.animation as anim

function csle2d(N, T, t0, dt, s, D, mu0, mu1)
    # initialize CSLE system
    z = zeros(ComplexF64,N,N)
    z2 = zeros(ComplexF64,N,N)
    dz = zeros(ComplexF64,N,N)
    s_sqrt_dt = s*sqrt(dt)
    Z = zeros(ComplexF64,T,N,N)
	# bifurcation parameter time course
    u = mu0.*ones(t0)
    v = [_ for _ in range(mu0,stop=mu1,length=T)]
    mu = vcat(u,v) # 1-dim time course of mu
    # iterate
    for t in range(1, stop=t0+T)
        (t%100 == 0) && print("    t = ", t, "\r")
        # CSLE equations
        z2 = z .* conj(z)
        dz = (mu[t] .+ 1im).*z - z.*z2 + D.*L(z)
        # stochastic integration
        z += (dz*dt + s_sqrt_dt .* randn(ComplexF64,N,N))
        (t > t0) && (Z[t-t0,:,:] = z)
    end
    println("\n")
    return Z
end

function animate_pyplot(fname, data)
    """
    Animate 3D array as .mp4 using PyPlot, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    vmin = minimum(data)
    vmax = maximum(data)
    # setup animation image
    println("[+] Animate")
    fig = figure(figsize=(6,6))
    axis("off")
    t = imshow(data[1,:,:], origin="lower", cmap=ColorMap("gray"),
               vmin=vmin, vmax=vmax)
    tight_layout()
    # frame generator
    function animate(i)
        (i%100 == 0) && print("    t = ", i, "/", nt, "\r")
        t.set_data(data[i+1,:,:])
    end
    # create animation
    ani = anim.FuncAnimation(fig, animate, frames=nt, interval=10)
    println("\n")
    # save animation
    ani[:save](fname, bitrate=-1,
               extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    show()
end

function animate_video(fname, data, downsample=10)
    """
    Animate 3D array as .mp4 using VideoIO, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    n1 = nt
    if (downsample > 0)
        println("[+] Downsampling ", downsample)
        n1 = Int(nt/downsample) # number of samples after downsampling
        data1 = zeros(Float64,n1,nx,ny)
        t0 = range(1,stop=nt)
        t1 = range(1,stop=nt,length=n1)
        for i in range(1,stop=nx)
            for j in range(1,stop=ny)
                y = data[:,i,j]
                f_ip = LinearInterpolation(t0, y)
                data1[:,i,j] = f_ip(t1)
            end
        end
    else
        data1 = copy(data)
    end
    # BW
    y = UInt8.(round.(255*(data .- minimum(data)) /
                          (maximum(data)-minimum(data))))
    # BW inverted
    #y = UInt8.(round.(255 .- 255*(data .- minimum(data)) /
    #                             (maximum(data)-minimum(data))))
    encoder_options = (color_range=2, crf=0, preset="medium")
    framerate=30
    T = size(data,1)
    println("[+] Animate")
    open_video_out(fname, y[1,end:-1:1,:], framerate=framerate,
                   encoder_options=encoder_options) do writer
        for i in range(2,stop=T,step=1)
            (i%100 == 0) && print("    i = ", i, "/", T, "\r")
            write(writer, y[i,end:-1:1,:])
        end
    end
    println("\n")
end

function L(x)
    # Laplace operator
    # periodic boundary conditions
    xU = circshift(x, [-1 0])
    xD = circshift(x, [1 0])
    xL = circshift(x, [0 -1])
    xR = circshift(x, [0 1])
    Lx = xU + xD + xL + xR - 4x
    # non-periodic boundary conditions
    Lx[1,:] .= 0.0
    Lx[end,:] .= 0.0
    Lx[:,1] .= 0.0
    Lx[:,end] .= 0.0
    return Lx
end

function main()
    println("Complex Stuart-Landau Equation (CSLE) lattice model\n")
    N = 128
    T = 2500
    t0 = 0
    dt = 0.05
    s = 0.05
    D = 1.0
    mu0 = -0.05
    mu1 = 0.5
    println("[+] Lattice size N: ", N)
    println("[+] Time steps T: ", T)
    println("[+] Warm-up steps t0: ", t0)
    println("[+] Integration time step dt: ", dt)
    println("[+] Noise intensity s: ", s)
    println("[+] Diffusion coefficient D: ", D)
    println("[+] CSLE parameter mu0: ", mu0, ", mu1: ", mu1)

    # run simulation
    Z = csle2d(N, T, t0, dt, s, D, mu0, mu1)
    data = angle.(Z)
    println("[+] Data dimensions: ", size(Z))

	# plot mean voltage
    m = mean(reshape(abs.(Z), (T,N*N)), dims=2)
    plot(m, "-k"); show()

	# save data
    mu0_str = rpad(mu0, 4, '0') # bif. parameter as 4-char string
    mu1_str = rpad(mu1, 4, '0') # bif. parameter as 4-char string
    s_str = rpad(s, 4, '0') # noise as 4-char string
    D_str = rpad(D, 4, '0') # diffusion coefficient as 4-char string
    fname1 = string("csle2d_mu0_", mu0_str, "_mu1_", mu1_str, "_s_", s_str,
					"_D_", D_str, ".npy")
    #npzwrite(fname1, data)
    #println("[+] Data saved as: ", fname1)

	fname2 = string("csle2d_mu0_", mu0_str, "_mu1_", mu1_str, "_s_", s_str,
					"_D_", D_str, ".mp4")
    animate_video(fname2, data, 4)
    println("[+] Data saved as: ", fname2)
end

main()


text:
module MyPackage

using Random

greet() = print("Hello World!")
hello(io::IO, name) = print(io, "Hello ", name)
hello(name) = hello(stdout, name)
randstr() = randstring()

end # module


text:
"""
Names of datasets and other objects used in simulation:

**Keys**
* `:osm` : an open street map object **(mandatory)**
* `:features` : an array of csv file names with informations about all the features (points representing crucial objects, like schools, recreation areas, shops locations, etc.) used in simulation **(mandatory)**
* `:flows` : csv file with informations about flows from each DA to another
* `:DAs` : csv file with centroids coordinates of every DA **(mandatory)**
* `:demo_stats` : csv file with demographic data of every DA **(mandatory)**
* `:googleapi_key` : a string containg unique *Google Distances API* key
"""
const filenames = Dict{Symbol,Union{String,Array{String,1}}}(:osm => "Winnipeg CMA.osm",
:flows =>"hw_flows.csv",
:DAs => "df_DA_centroids.csv",
:googleapi_key => "googleapi.key"
)

"""
A dictionary of dictionaries structuring the way how the demographic data will influence the demographic profile of agents in simulation. It is based on the mandatory `:demo_stats` data frame which is containing an informations about socio-demographic structure of each DA.

"""
const demographic_categories = Dict(
    :gender => Dict(
    :male => "Male",
    :female => "Female"
    ),

    :age => Dict(
    :age20to29 => (20:29), #"Household Population by Age - 20 To 29",
    :age30to39 => (30:39), #"Household Population by Age - 30 To 39",
    :age40to49 => (40:49), #"Household Population by Age - 40 To 49",
    :age50to59 => (50:59), #"Household Population by Age - 50 To 59",
    :age60to69 => (60:69), #"Household Population by Age - 60 To 69",
    ),

    :marital_status => Dict(
    :married  => "Married",
    :unmarried => "Unmarried"
    ),

    :household_income => Dict(
    :income_19999_or_less => (0:19999),        #"Income - 19,999 Or Less",
    :income_20000_to_39999 => (20000:39999),   #"Income - 20,000 To 39,999",
    :income_40000_to_59999 => (40000:49999),   #"Income - 40,000 To 59,999",
    :income_60000_to_79999 => (60000:79999),   #"Income - 60,000 To 79,999",
    :income_80000_to_99999 => (80000:99999),   #"Income - 80,000 To 99,999",
    :income_80000_to_124999 => (100000:124999), #"Income - 100,000 To 124,999",
    :income_125000_to_149999 => (125000:149999), #"Income - 125,000 To 149,999",
    :income_150000_to_174999 => (150000:174999), #"Income - 150,000 To 174,999",
    :income_175000_to_199999 => (175000:199999), #"Income - 175,000 To 199,999",
    :income_200000_to_249999 => (200000:249999), #"Income - 200,000 To 249,999",
    :income_250000_or_over => (250000:1000000) #"Income - 250,000 Or Over",
    ),

    :no_of_children => Dict(
    :nochildren => 0,
    :onechild  => 1,
    :twochildren  => 2,
    :threechildrem  => 3,
    ),

)


text:
B_data = [ -2.5           0.5           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0          -2.0           1.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
1.0          -3.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.0          -2.0           1.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0          -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.5           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0          -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           2.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.5          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           2.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.5           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           2.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.5           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   2.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.5          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   2.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.5           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   2.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.5        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   2.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0        0.5        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   2.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.5       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.5           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           2.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.5          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           2.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.5           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0   0.0   0.0   0.0           0.0           0.0           2.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.5           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   2.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.5          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   2.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.5           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0  -0.0   0.0           0.0           0.0           0.0   0.0   0.0   2.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.5        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   2.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0        0.5        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -2.5   0.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   2.0   0.0           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.5       0.0        0.0
0.0           0.0           1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -5.0           5.0          -2.0           0.0   0.0   0.0   0.0   0.0   0.0   1.0           0.0          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.5          -2.0           0.5           0.0   0.0   0.0   0.0   0.0   0.0   0.0           1.0          0.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.5          -2.0          -3.5           4.5  -3.0   1.5   0.0   0.0   0.0   0.0           0.0          1.0           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0          -5.0   5.0  -2.0   0.0   0.0   0.0   0.0           0.0          0.0           1.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0          -0.5  -2.0   0.5   0.0   0.0   0.0   0.0           0.0          0.0           0.0           1.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.5  -2.0  -3.5   4.5  -3.0   1.5   0.0           0.0          0.0           0.0           0.0          1.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0  -5.0   5.0  -2.0   0.0           0.0          0.0           0.0           0.0          0.0           1.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0  -0.5  -2.0   0.5   0.0           0.0          0.0           0.0           0.0          0.0           0.0        1.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.5  -2.0  -3.5   0.0           0.0          0.0           0.0           0.0          0.0           0.211839  -0.423678   1.63552   0.436525   2.1398
3.26466       0.683008      0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.8285       -0.552334      0.276167      0.0   0.0   0.0   0.0   0.0   0.0  -5.25         -3.0          0.75          0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
1.35919e-16   2.84359e-17   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   3.44933e-17   2.0           1.14978e-17   0.0   0.0   0.0   0.0   0.0   0.0   0.75         -3.0         -0.75          0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
-9.06125e-16  -1.89573e-16   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.29955e-16   1.53303e-16   2.0           0.0   0.0   0.0   0.0   0.0   0.0  -3.0           7.5         -7.5           0.0           0.0          0.0           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           2.0   0.0   0.0   0.0   0.0   0.0   2.25         -4.5          6.75         -5.25         -3.0          0.75          0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   2.0   0.0   0.0   0.0   0.0   9.36751e-17  -1.8735e-16   2.81025e-16   0.75         -3.0         -0.75          0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   2.0   0.0   0.0   0.0  -6.245e-16     1.249e-15   -1.8735e-15   -3.0           7.5         -7.5           0.0        0.0        0.0       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   2.0   0.0   0.0   0.0           0.0          0.0           2.25         -4.5          6.75         -5.25      -3.0        0.75      0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   2.0   0.0   0.0           0.0          0.0           9.36751e-17  -1.8735e-16   2.81025e-16   0.75      -3.0       -0.75      0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   2.0   0.0           0.0          0.0          -6.245e-16     1.249e-15   -1.8735e-15   -3.0        7.5       -7.5       0.0        0.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.25      -0.5        0.75     -2.5        2.0
0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0           0.0           0.0   0.0   0.0   0.0   0.0   0.0   0.0           0.0          0.0           0.0           0.0          0.0           0.5       -1.0        1.5       1.0       -2.0]


text:
module MTKHelpers

using ModelingToolkit, DifferentialEquations
using StaticArrays, LabelledArrays
using Requires: @require
using NamedArrays
using ComponentArrays
#using Infiltrator

export symbol, symbols_state, symbols_par, strip_namespace, embed_system

export AbstractProblemParSetter, ProblemParSetter_sym, ProblemParSetter,
    count_state, count_par, count_paropt,
    axis_paropt, axis_par, axis_state,
    symbols_paropt, symbols_state, symbols_par,
    update_statepar, get_paropt, get_paropt_labeled, label_paropt, label_par, label_state,
    name_paropt, name_par, name_state

export
    _get_index_axis, _set_index_axis!, attach_axis, _update_cv, _labels

export smoothstep

export series_sol!

export getlast

# extending
import Base: merge, getindex
import ComponentArrays: getdata

function __init__()
    @require CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" include("requires_cairommakie.jl")
  end

include("util.jl")
include("util_componentarrays.jl")
include("abstractproblemparsetter.jl")
include("problemparsetter_sym.jl")
include("problemparsetter.jl")
include("smoothstep.jl")
include("solution.jl")

end


text:
"""
This is a docstring! If you do ?added_function in the REPL
or in your IDE, this string will be the printed help string.

```julia
x + y
```

Inline code can be added as well, and math (LaTeX) can be
done in a

```math
\Vert x - y \Vert
```

math block. Docstrings can be applied to more than just
functions: types, globals, etc. can all have docstrings.
Documenter.jl can use these docstrings in your documentation
as well.
"""
added_function(x,y) = x+y


text:
module GrayCoding

using LinearAlgebra
using Gadfly

export πmatrix

# Write your package code here.
"""
Generate Encoding and Decoding matrices for Gray Codes of alphabet.
```julia-repl
julia> G,B,g,b=GrayMatrix(4, 2);
julia> G
    4×4 Matrix{Int64}:
    1  0  0  0
    1  1  0  0
    1  1  1  0
    1  1  1  1
    julia> B
    4×4 Matrix{Int64}:
    1  0  0  0
    1  1  0  0
    0  1  1  0
    0  0  1  1
    julia> g
    4×16 Matrix{Int64}:
    0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
    0  0  0  0  1  1  1  1  0  0  0  0  1  1  1  1
    0  0  1  1  0  0  1  1  0  0  1  1  0  0  1  1
    0  1  0  1  0  1  0  1  0  1  0  1  0  1  0  1
    julia> b
    4×16 Matrix{Int64}:
    0  0  0  0  0  0  0  0  1  1  1  1  1  1  1  1
    0  0  0  0  1  1  1  1  1  1  1  1  0  0  0  0
    0  0  1  1  1  1  0  0  0  0  1  1  1  1  0  0
    0  1  1  0  0  1  1  0  0  1  1  0  0  1  1  0
```
"""
function GrayMatrix(n::Int64=3, q::Int64=2)
    M = q^n
    G = Array(Bidiagonal(ones(1, n)[:], (q - 1) * ones(1, n - 1)[:], :L))
    B = mod.(inv(G), q)
    G = convert.(Int64, G)
    B = convert.(Int64, B)
    x = 0:M-1
    a = digits.(x, base = q, pad = n)  # left-lsb
    # a=digits.(x,base=2,pad=num_bits)|> reverse  # left-msb
    b = hcat(a...)
    b = reverse(b, dims = 1)
    g = Int.(mod.(G * b, q))
    return B, G, b, g
end

"""
Plots a matrix into a 2D with labels. Optional arguments including colors
```julia-repl
julia> using Random;
julia> A= rand(0:9,10,10);
julia> matrixplot(A)
```
"""
function matrixplot(A;kwargs...)
    a,b=size(A)
    X = transpose(repeat(1:b, 1, a))[:]
    Y = repeat(a:-1:1, b)[:]
    Gadfly.plot(x = X, y = Y, Geom.rectbin(),  color = A,alpha=[0.5], Coord.cartesian(fixed = true), Theme(bar_spacing = 0.1cm), Geom.label(position = :centered), label = string.(A)[:], Theme(key_position = :none, grid_line_width = 0pt, panel_stroke = nothing), Guide.xticks(ticks = nothing, label = false), Guide.yticks(ticks = nothing, label = false), Guide.xlabel(nothing), Guide.ylabel(nothing);kwargs...)
end

"""
Plot the DNA codon matrix
```julia-repl
julia> dnamatrix()
```
"""
function dnamatrix()
	U = [(x,y) for x ∈ 0:7, y ∈ 0:7]
	B,G,b,g=GrayMatrix(3,2)
	V = [string(g[:,y+1]...)string(g[:,x+1]...) for x ∈ 0:7, y ∈ 0:7]
	revV = [string(g[:,y+1]...)string(g[:,x+1]...) for x ∈ 0:7, y ∈ 0:7]
	Vx = [string(g[:,y+1]...) for x ∈ 0:7, y ∈ 0:7]
	Vy = [string(g[:,x+1]...) for x ∈ 0:7, y ∈ 0:7]
	VM=[(Vx[i,j][1]Vy[i,j][1]Vx[i,j][2]Vy[i,j][2]Vx[i,j][3]Vy[i,j][3]) for i ∈ 1:8, j ∈ 1:8]
	DNA=parse.(Int,VM,base=2)
	# DM0=[replace(Vx[i,j][1]Vy[i,j][1],"00"=>"C","01"=>"A","10"=>"U","11"=>"G")replace(Vx[i,j][2]Vy[i,j][2],"00"=>"C","01"=>"A","10"=>"U","11"=>"G")replace(Vx[i,j][3]Vy[i,j][3],"00"=>"C","01"=>"A","10"=>"U","11"=>"G") for i ∈ 1:8, j ∈ 1:8]
	DM=[replace(Vx[i,j][1]Vy[i,j][1],"00"=>"C","01"=>"A","10"=>"U","11"=>"G")replace(Vx[i,j][2]Vy[i,j][2],"00"=>"C","01"=>"A","10"=>"U","11"=>"G")replace(Vx[i,j][3]Vy[i,j][3],"00"=>"C","01"=>"A","10"=>"U","11"=>"G") for j ∈ 1:8, i ∈ 1:8]
    AA=copy(DM) # Amino Acid
	replace!(AA,"CUU"=>"Leucine","CUC"=>"Leucine","CUA"=>"Leucine","CUG"=>"Leucine","UUA"=>"Leucine","UUG"=>"Leucine","UUU"=>"Phenylalanine","UUC"=>"Phenylalanine","AUC"=>"Isoleucine","AUA"=>"Isoleucine","AUU"=>"Isoleucine","AUG"=>"Methionine","GUA"=>"Valine","GUC"=>"Valine","GUU"=>"Valine","GUG"=>"START","UCA"=>"Serine","UCC"=>"Serine","UCU"=>"Serine","UCG"=>"Serine","CCC"=>"Proline","CCA"=>"Proline","CCU"=>"Proline","CCG"=>"Proline","ACU"=>"Threonine","ACA"=>"Threonine","ACC"=>"Threonine","ACG"=>"Threonine","GCC"=>"Alanine","GCU"=>"Alanine","GCA"=>"Alanine","GCG"=>"Alanine","GGU"=>"Glycine","GGA"=>"Glycine","GGC"=>"Glycine","GGG"=>"Glycine","CGA"=>"Arginine","CGC"=>"Arginine","CGU"=>"Arginine","CGG"=>"Arginine","GAU"=>"Aspartic acid","GAC"=>"Aspartic acid","GAA"=>"Glutamic acid","GAG"=>"Glutamic acid","AAU"=>"Asparagine","AAC"=>"Asparagine","UGU"=>"Cysteine","UGC"=>"Cysteine","UGG"=>"Tryptophan","CAA"=>"Glutamine","CAG"=>"Glutamine","UAA"=>"STOP","UAG"=>"STOP","UAU"=>"Tyrosine","UAC"=>"Tyrosine","AAA"=>"Lysine","AAG"=>"Lysine","CAC"=>"Histidine","CAU"=>"Histidine","AGG"=>"Arginine","AGA"=>"Arginine","AGU"=>"Serine","AGC"=>"Serine","UGA"=>"STOP" )
	return DM,VM,AA,Vx,Vy,DNA
end

"""
Decimal to binary conversion
```julia-repl
julia> dec2bin(10,7)
```
"""
function dec2bin(x,n)
	a=digits.(x,base=2,pad=n)  # left-lsb
# a=digits.(x,base=2,pad=num_bits)|> reverse  # left-msb
b=hcat(a...);
	return b
end

"""
Binary to decimal number conversion. Input can be
-   binary strings,
-   binary digits or
-   a vector of binary string or digits
```julia-repl
julia> bin2dec([011,"0111",0111])
```
"""
bin2dec(u) =[parse(Int,string(u[:,ii]...),base=2) for ii in 1:size(u,2)]
bin2dec(x::AbstractVector)  = parse(Int,string(x...),base=2)
bin2dec(x::Int64)  = parse(Int,string(x),base=2)
bin2dec(x::AbstractString) = parse(Int,x,base=2)

"""
Pulse amplitude modulation (PAM) mapping. This is a type of digital modulation mapping used in Communication systems.

"""
function pam_encode(x,M)
	# M --> M-QAM
	n=Int(ceil(log2(M)/2))
	B,G,b,g=find_gray_matrix(n)
	u=digits(x,base=2,pad=n) |> reverse
	v=Int.(mod.(G*u,2))
	w=bin2dec(v)
	y=-sqrt(M)+1+2*w
	return y
end

"""
Generate Gray vectors
"""
function gen_gray(m)
	x=[0 1]
	for i in 2:m
		rev_x=reverse(x)
		x=[x rev_x.+2^(i-1)]
	end
	return x[:]
end

"""
Decimal to binary conversion
"""
function dec2bin(x,n)
	a=digits.(x,base=2,pad=n)  # left-lsb
# a=digits.(x,base=2,pad=num_bits)|> reverse  # left-msb
    b=hcat(a...);
	return b
end

"""
Recursive function to illustrate the reflection+shift property of Gray mapping.
### Arguents
* n - The iteration number `n ≥ 0`
* C - The decimal sequence of the gray mapped bits
* R - The reflected sequence (without shifting)
```julia-repl
julia> C,R = gray_recursion(4)
```
"""
function gray_recursion(n::Int)
    C= n < 1 ? [0] : vcat(gray_recursion(n - 1)[1:end], gray_recursion(n - 1)[end:-1:1] .+ Int(exp2(n-1)) )
	return C
end

"""
Reflected code.
```julia-repl
julia>reflect_code(3)
[0,1,3,2,2,3,1,0]
```
"""
function reflect_code(n)
    n<1 ? [0] : vcat(gray_recursion(n-1),reverse(gray_recursion(n-1)))
end

"""
Recursive construction of binary Gray code digits.

Gray code ``g[n]`` can be recursively constructed as follows.
Start with `` g[1] = (0,1) = (g_{1},g_{2})``

```math
g[n+1] = 0g_{1},...,0g_{N−1},0g_{N},1g_{N},1g_{N−1},...,1g_{1}.
```

### Examples
```julia-repl
julia> gray(3)
3×8 Matrix{Int64}:
 0  0  0  0  1  1  1  1
 0  0  1  1  0  1  1  0
 0  1  1  0  1  1  0  0
```
"""
function gray(n)
    n < 2 ? [0 1] : hcat(vcat(0,gray(n-1)),vcat(1,reverse(gray(n-1))))
end

"""
Find the Givens embedding ``{}^{i}G_{j,k}(A)``
"""
function G(A,i,j,k)
	n = size(A,2)
	G = Matrix(1.0I(n).+0im)
	α = A[k,i]
	β = A[j,i]

	Γ0 = [ α' β'
	      -β  α ]
	N = norm(Γ0[1,:],2);
	Γ = Γ0./N

    # Embed the Givens matrix Γ in G
	G[k,k] = Γ[1,1]
	G[k,j] = Γ[1,2]
	G[j,k] = Γ[2,1]
	G[j,j] = Γ[2,2]

	return G,G*A
end

"""
Find the matrix which has a Givens matrix embedding ``{}^{i}G_{j,k}(A)``.

For a given matrix ``A``, a generic rotation matrix ``{}^{i}G_{j,k}(A)`` is generated. The matrix ``{}^{i}G_{j,k}(A)`` is such that it helps to selectively nullifies an element of matrix ``V={}^{i}G_{j,k}(A) A``. That is, ``V_{j,i}=0`` where ``V={}^{i}G_{j,k}(A) A``.  The fllowing Givens rotation matrix ``{}^{i}\\Gamma_{j,k}``

```math
\\begin{aligned}
{}^{i}\\Gamma_{j,k} &= \\begin{pmatrix} {}^{i}g_{k,k} & {}^{i}g_{k,j} \\\\ {}^{i}g_{j,k} & {}^{i}g_{j,j}\\end{pmatrix} \\\\
&=\\frac{1}{\\sqrt{\\lvert a_{j,i} \\rvert^{2}+ \\lvert a_{k,i} \\rvert^{2}}}\\begin{pmatrix} a_{k,i}^{*} & a_{j,i}^{*} \\\\ -a_{j,i} & {}^{i}a_{k,i}\\end{pmatrix}.
\\end{aligned}
```

is embedded in an identity matrix ``I(n)`` to produce,
```math
{}^{i}G_{j,k}(A) = \\begin{pmatrix}
1 & 0 & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & 0 \\\\
0 & 1 & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\
\\vdots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\
0 & 0 & \\ddots & {\\color{red} {}^{i}g_{k,k}} & \\ddots &  {\\color{red} {}^{i}g_{k,j}} & \\ddots & \\ddots & \\vdots \\\\
\\vdots & \\ddots & \\ddots & \\ddots & \\ddots \\ddots & & \\ddots & \\ddots & \\vdots \\\\
0 & 0 & \\ddots & {\\color{red}  {}^{i}g_{j,k}} & \\ddots &  {\\color{red}
 {}^{i}g_{j,j}} & \\ddots & \\ddots & \\vdots \\\\
\\vdots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\
0 & 0 & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & 1 & 0  \\\\
0 & 0 & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & 1
\\end{pmatrix}.
```
Essentially, ``{}^{i}G_{j,k}(A)`` is a modified identity matrix such that four non trivial elements are taken from the givens rotation matrix ``{}^{i}\\Gamma_{j,k}``.

"""
function GivensG(A,i,j,k)
	n = size(A,2)
	G = Matrix(1.0I(n).+0im)
	α = A[k,i]
	β = A[j,i]

	Γ0 = [ α' β'
	      -β  α ]
	N = norm(Γ0[1,:],2);
	Γ = Γ0./N

    # Embed the Givens matrix Γ in G
	G[k,k] = Γ[1,1]
	G[k,j] = Γ[1,2]
	G[j,k] = Γ[2,1]
	G[j,j] = Γ[2,2]

	return G,G*A
end

"""
Given unitary matrux ``U(n) \\in S(2^n)``, it uses a a repeated Givens rotations to get a to level matrix as follows.

```math
\\prod_{j=1}^{n-2}{ {}^{1}G_{n-j,n-j-1} U(n)}  = \\begin{pmatrix} 1 & 0 \\\\ 0 & U(n-1)\\end{pmatrix}
```

### Parameters
* U -- Input. Unitary matrix of size ``2^n``
* V -- Output unitary matrix in two level form `[1 0;0 U']` form where `U'` is a unitary matrix of size ``2^{n-1}``.
* GG -- The ``n-2`` sequence of Given matrices (in augmented form) ``[{}^{1}G_{n}\\lvert{}^{1}G_{n-1}\\lvert\\ldots\\lvert{}^{1}G_{2}]``

### Examples
```julia-repl
julia> level2unitary(U)
```
"""
function level2unitary(U)
	n=size(U,2)
	V=I(n)*U;
	GG=Matrix(I(n))
	for i=1:1
		for j=0:n-1-i
			G0,V=GivensG(V,i,n-j,n-j-1)
			GG=[GG;G0]
		end
	end
	return GG[n+1:end,:],V
end

"""
Decomposition of aribtrary unitary matrix ``U(n) \\in S(2^n)``, as a cascade of two level Givens matrices.

Using the following property,
```math
\\prod_{j=0}^{n-1}{ {}^{1}G_{n-j,n-j-1} U(n)}  = \\begin{pmatrix} 1 & 0 \\\\ 0 & U(n-1)\\end{pmatrix}
```

```math
\\prod_{i=1}^{n-1} \\prod_{j=0}^{n-1-i}{ {}^{i}G_{n-j,n-j-1} U(n)}  = \\begin{pmatrix} 1 & \\ldots & 0 \\\\ \\vdots & \\ddots & \\vdots \\\\ 0 & \\ldots & \\det(U)\\end{pmatrix}.
```
There are ``2^{n-1}(2^{n}-1)`` number of unitary two-level matrices (each of them formed by embedding a Givens rotaton matrix into indentity matrix). Note that ``\\sum_{i=1}^{2^{n}}{N-i}=2^{n-1}(2^{n}-1)``.

### Parameters
* U -- Input. Unitary matrix of size ``2^n``
* V -- Output unitary matrix in two level form `[1 0;0 U']` form where `U'` is a unitary matrix of size ``2^{n-1}``.
* Gm -- The ``(n-2)(n-1)`` sequence of Given matrices (in augmented form) ``[{}^{1}G_{n}\\lvert{}^{1}G_{n-1}\\lvert\\ldots\\lvert{}^{1}G_{2}\\lvert {}^{2}G_{n}\\lvert{}^{2}G_{n-1}\\lvert\\ldots\\lvert{}^{2}G_{3}\\lvert\\ldots\\lvert{}^{n-2}G_{n}\\lvert{}^{n-2}G_{n-1}\\lvert{}^{n-1}G_{n}]``
* Gs -- The ``(n-2)(n-1)`` sequence of Given matrices (in augmented form) left to right ``[ {}^{n-1}G_{n} \\lvert {}^{n-2}G_{n-1} \\lvert {}^{n-2}G_{n} \\lvert \\ldots \\lvert {}^{2}G_{3} \\lvert \\ldots \\lvert {}^{2}G_{n-1} \\lvert {}^{2}G_{n}  \\lvert{}^{1}G_{2} \\lvert \\ldots \\lvert {}^{1}G_{n-1} \\lvert {}^{1}G_{n} ]``

### Examples
```julia-repl
julia> using LinearAlgebra
julia> N=4;A=rand(N,N)+im*rand(N,N);S=svd(A);U=S.U
julia> Gn(U)
```
"""
function Gn(A)
	n=size(A,2)
	V=I(n)*A;
	Gm=Matrix(I(n))
	Gs=Matrix(I(n))
	for i=1:n-1
		for j = 0:n-1-i
			G1,V = GivensG(V,i,n-j,n-j-1)
			Gm   = [Gm ; G1]
			Gs   = [G1 ; Gs]
		end
	end
	return V,Gm[n+1:end,:],Gs[1:end-n,:]

end

"""
For a given unitary matrix ``U``, it finds the cascaded rotation matrix ``C`` such that ``C\\times U =I``, except for the diagonal element of the ``n``th element which is ``\\det(U)``. The matrix ``C`` is obtained by the cascade of several (i.e., ``2^{n-1}(2^{n}-1)``) two level Givens rotation matrices, namely,

```math
C  = \\prod_{i=1}^{n-1} \\prod_{j=0}^{n-1-i}{ {}^{i}G_{n-j,n-j-1}}
```

```math
\\prod_{i=1}^{n-1} \\prod_{j=0}^{n-1-i}{ {}^{i}G_{n-j,n-j-1} U(n)}  = \\begin{pmatrix} 1 & \\ldots & 0 \\\\ \\vdots & \\ddots & \\vdots \\\\ 0 & \\ldots & \\det(U)\\end{pmatrix}.
```

### Parameters
* U -- Input. Unitary matrix of size ``n``
* Ic --  Identity matix with the exception that the last diagonal entry (``n``th diagonal element) which is ``\\det(U)``.
* C -- The cascaded rotation matrix ``{}^{n-1}G_{n} \\times {}^{n-2}G_{n-1}  \\times {}^{n-2}G_{n}  \\times \\ldots  \\times {}^{2}G_{3}  \\times \\ldots  \\times {}^{2}G_{n-1}  \\times {}^{2}G_{n}   \\times {}^{1}G_{2}  \\times \\ldots  \\times {}^{1}G_{n-1}  \\times {}^{1}G_{n} ``

### Examples
```julia-repl
julia> using LinearAlgebra
julia> N=4;A=rand(N,N)+im*rand(N,N);S=svd(A);U=S.U
4×4 Matrix{ComplexF64}:
  -0.4-0.06im   0.23-0.73im  -0.03-0.14im   0.39-0.27im
 -0.61-0.32im   0.07+0.06im   0.32+0.02im  -0.64-0.08im
 -0.38-0.33im  -0.48+0.38im  -0.07-0.4im    0.46+0.07im
  -0.2-0.26im  -0.09-0.14im   -0.7+0.47im  -0.09+0.38im
julia> Gc,Ic=Gcascaded(U)
4×4 Matrix{ComplexF64}:
 1.0+0.0im   0.0+0.0im  -0.0+0.0im   0.0-0.0im
 0.0+0.0im   1.0-0.0im   0.0+0.0im  -0.0-0.0im
 0.0-0.0im  -0.0+0.0im   1.0+0.0im   0.0+0.0im
 0.0-0.0im  -0.0+0.0im   0.0-0.0im   0.4-0.9im
julia> det(A)
0.4166084175706718 - 0.9090860390575042im
```
"""
function Gcascaded(U)
	Gv=Gn(U)[3]
	Gc=Matrix(I(4))
	m=Int(size(Gv,1)/4)
	for q in 0:m-1
		Gc=Gc*Gv[1+q*4:4*(q+1),:]
	end
	Ic=Gc*U
	return Gc,Ic
end

"""
For ``n`` bits, with the corresponding decimal sequence ``x=0,1,2,\\ldots,2^{n-1}``, find the gray ordering sequence using the bitwise `XOR` logic. Namely.
``gπ= x ⊻ \\lfloor x/2 \\rfloor = x \\oplus \\lfloor x/2 \\rfloor`` where ``⊻`` is the bitwise `XOR` operation.

```julia-repl
julia> gray_ordering(3)
[0,1,3,2,6,7,5,4]
```
"""
function gray_ordering(n)
	N=2^n
	x=Int.(vcat(range(start=0,stop=N-1,step=1)))
	y = Int.(floor.(x/2))
	gπ=[x[i] ⊻ y[i] for i=1:N]
	return gπ
end

""" Converts a gray mapped decimal number to simple binary mapped decimal number. This routine follows the inverse mapping of the function `gray_ordering()`
### Parameters
* ``n`` The number of digits in the equivalent binary mapping
* ``g`` The gray mapped decimal number
* ``d`` The simple binary mapped decimal number

```julia-repl
julia> G2D.(0:7,3)
[0,1,3,2,7,6,4,5]

```

"""
function G2D(g::Int,n)
	xx=gray_ordering(n) .+1
	a=collect(1:2^n)
	p=invperm(xx)
	d=p[g+1]
	return d-1
end

"""
Quantum circuit decomposition. For nan arbitrary unitary matix for ``n`` qubits.
Arbitrary quantum circuit abstracted by unitary matrix ``U``  decomposed by ``2^{n-1}2^{n}``
unitary two-level matrices, each of which corresponds ``{}^{i}\\Gamma_{j,k}``.
The program produce the ``\\Gamma`` matrix and the coefficients ``i,j,k``. The quantum
circuit of this decomposition can be visualized as a cascade (from left to right) of this matrix.
"""
function sequenceΓ(k::Int64)
 	n=Int(2^k)
	N=Int(n*(n-1)/2)
	idxI=-100
	idxJ=-100
	m=0
	for i=1:n-1
		for j = n:-1:(i+1)
			idxI=[idxI n-i]
			idxJ=[idxJ j]

		end
	end
	dI,dJ=Int.(idxI[2:end] .+0),Int.(idxJ[2:end] .+0)
 	gV=gray_ordering(k)

	gJ2 = gV[dJ .- 1] .+1

	gI=gV[-dI.+0 .+ n] .+1

	gJ= gV[dJ] .+1

	ii=reverse(gI)
	jj=reverse(gJ)
	kk=reverse(gJ2)
	Γ = hcat(ii,jj,kk)'
	return ii,jj,kk,Γ

end

"""
The unitary matrix ``U`` corresponding to of the 3 quibit Tiffoli quantum gate

```math
U=\\begin{pmatrix}
1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\\\
0 & 0 & 0 & 0 & 0 & 0 & {\\color{red}0} & {\\color{red}1} \\\\
0 & 0 & 0 & 0 & 0 & 0 & {\\color{red}1} & {\\color{red}0} \\\\
\\end{pmatrix}
```
"""
function toffoli_matrix()
	U=Matrix(I(8))
	U[8,8]=U[7,7]=0
	U[8,7]=U[7,8]=1
	return U
end

"""
Given a matrix ``A`` (usually a unitary matrix) and indices ``(i,j,k)``, find the ``\\Gamma`` Givens rotation matrix and the Givens matrix embedding ``{}^{i}G_{j,k}(A)`` matrix.

For a given matrix ``A``, a generic rotation matrix ``{}^{i}G_{j,k}(A)`` is generated. The matrix ``{}^{i}G_{j,k}(A)`` is such that it helps to selectively nullifies an element of matrix ``V={}^{i}G_{j,k}(A) A``. That is, ``V_{j,i}=0`` where ``V={}^{i}G_{j,k}(A) A``.  The fllowing Givens rotation matrix ``{}^{i}\\Gamma_{j,k}``

```math
\\begin{aligned}
{}^{i}\\Gamma_{j,k} &= \\begin{pmatrix} {}^{i}g_{k,k} & {}^{i}g_{k,j} \\\\ {}^{i}g_{j,k} & {}^{i}g_{j,j}\\end{pmatrix} \\\\
&=\\frac{1}{\\sqrt{\\lvert a_{j,i} \\rvert^{2}+ \\lvert a_{k,i} \\rvert^{2}}}\\begin{pmatrix} a_{k,i}^{*} & a_{j,i}^{*} \\\\ -a_{j,i} & {}^{i}a_{k,i}\\end{pmatrix}.
\\end{aligned}
```

is embedded in an identity matrix ``I(n)`` to produce,
```math
{}^{i}G_{j,k}(A) = \\begin{pmatrix}
1 & 0 & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & 0 \\\\
0 & 1 & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\
\\vdots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\
0 & 0 & \\ddots & {\\color{red} {}^{i}g_{k,k}} & \\ddots &  {\\color{red} {}^{i}g_{k,j}} & \\ddots & \\ddots & \\vdots \\\\
\\vdots & \\ddots & \\ddots & \\ddots & \\ddots \\ddots & & \\ddots & \\ddots & \\vdots \\\\
0 & 0 & \\ddots & {\\color{red}  {}^{i}g_{j,k}} & \\ddots &  {\\color{red}
 {}^{i}g_{j,j}} & \\ddots & \\ddots & \\vdots \\\\
\\vdots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & \\vdots \\\\
0 & 0 & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & 1 & 0  \\\\
0 & 0 & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & \\ldots & 1
\\end{pmatrix}.
```
Essentially, ``{}^{i}G_{j,k}(A)`` is a modified identity matrix such that four non trivial elements are taken from the givens rotation matrix ``{}^{i}\\Gamma_{j,k}``.
```julia-repl
julia> using LinearAlgebra
julia> using GrayCoding
julia>n=2;N=2^n;A=rand(N,N)+im*rand(N,N);S=svd(A);U=S.U
4×4 Matrix{ComplexF64}:
 -0.365903-0.405021im   0.442293-0.0769938im  …   0.115307-0.288609im
 -0.285173-0.35669im   -0.671764+0.0698449im     -0.384583+0.295428im
 -0.196831-0.611652im  -0.154487+0.0160399im      0.379159-0.121825im
 -0.177839-0.221435im   0.536228-0.175044im      -0.338822+0.62835im
julia> i,j,k=1,2,4
julia> Γ,G,GA=quantumΓ(U,i,j,k);
julia> round.(quantumΓ(S.U,1,2,4)[1],digits=1)
2×2 Matrix{ComplexF64}:
 -0.3+0.4im  -0.5+0.7im
  0.5+0.7im  -0.3-0.4im
julia> round.(quantumΓ(S.U,1,2,4)[2],digits=1)
4×4 Matrix{ComplexF64}:
 1.0+0.0im   0.0+0.0im  0.0+0.0im   0.0+0.0im
 0.0+0.0im  -0.3-0.4im  0.0+0.0im   0.5+0.7im
 0.0+0.0im   0.0+0.0im  1.0+0.0im   0.0+0.0im
 0.0+0.0im  -0.5+0.7im  0.0+0.0im  -0.3+0.4im
julia> round.(quantumΓ(S.U,1,2,4)[3],digits=1)
4×4 Matrix{ComplexF64}:
 -0.4-0.4im   0.4-0.1im   0.6+0.1im   0.1-0.3im
  0000000     0.7+0.5im  -0.2-0.4im  -0.3+0.2im
 -0.2-0.6im  -0.2+0.0im  -0.4-0.5im   0.4-0.1im
  0.5+0.0im   0.2-0.2im  -0.2-0.1im  -0.1-0.8im
```

"""
function quantumΓ(A,i,j,k)
	n = size(A,2)
	G = Matrix(1.0I(n).+0im)
	α = A[k,i]
	β = A[j,i]

	Γ0 = [ α' β'
	      -β  α ]
	N = norm(Γ0[1,:],2);
	Γ = Γ0./N

    # Embed the Givens matrix Γ in G
	G[k,k] = Γ[1,1]
	G[k,j] = Γ[1,2]
	G[j,k] = Γ[2,1]
	G[j,j] = Γ[2,2]

	return Γ,G,G*A
end

"""
Permutation matrix ``P`` which when operated on the right will swap columns ``(k,n)``. When operated n the left the rows ``(k,n)`` will get swapped. Permutation matrices are derived from identity matrix by swapping one column (also results in one row swap).
##### Arguments
* ``n`` -- order of the matrix (square matrix of size ``n \\times n``)
* ``k,m`` -- column indices ``1 \\ge (k,m) \\ge n``. When ``k,m`` is a vector, multiple rows/columns gets swapped when operated on left/right.
##### Examples
```julia-repl
julia> πmatrix(6,[2,3],[5,6])
6×6 Matrix{Float64}:
 1.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  1.0
 0.0  0.0  0.0  1.0  0.0  0.0
 0.0  1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0  0.0
julia> A= rand(4,4)
4×4 Matrix{Float64}:
 0.335342   0.880332  0.747449  0.886538
 0.0326128  0.648394  0.87131   0.386876
 0.876305   0.215459  0.247227  0.100161
 0.707862   0.992129  0.411607  0.995662
julia> AA*πmatrix(4,2,4)
4×4 Matrix{Float64}:
 0.335342   0.886538  0.747449  0.880332
 0.0326128  0.386876  0.87131   0.648394
 0.876305   0.100161  0.247227  0.215459
 0.707862   0.995662  0.411607  0.992129

julia> πmatrix(4,1,3)*AA
4×4 Matrix{Float64}:
 0.876305   0.215459  0.247227  0.100161
 0.0326128  0.648394  0.87131   0.386876
 0.335342   0.880332  0.747449  0.886538
 0.707862   0.992129  0.411607  0.995662
```

#### Test
```julia
	begin
		R8=rand(8,8)
		Q8=πmatrix(8,2,7)*πmatrix(8,5,8)*R8*πmatrix(8,2,7)*πmatrix(8,5,8)
		R8[[2,5],[2,5]]
		Q8[[7,8],[7,8]]
	end
```

"""
function πmatrix(n=5,k=2,m=n-1)
	M=diagm(0=>ones(n))	# I(m)
	p=collect(1:n)
	p[k]=m
	p[m]=k
	P=M[:,p]
	return P
end

"""
Given an integer ``x ∈ \\mathbb{N}`` (i.e., natural number in decimal form), it performs a tensor product of the individual bits. Effectively, each of the bits in the binary representation expressed in a qubit vector gets tensor multipled.
The qubit to vector mapping are as follows:
``0 \\equiv | 0 \\rangle = \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix}`` and ``1 \\equiv | 1 \\rangle = \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix}``

##### Example
``x=13`` in binary representation is given by `` 1101 ``. Thus ``|13 \\rangle \\equiv |1101\\rangle=|1\\rangle \\otimes |1\\rangle \\otimes |0\\rangle \\otimes |1\\rangle``.
```math
\\begin{aligned}
|13 \\rangle &\\equiv |1101\\rangle \\\\
&=|1\\rangle \\otimes |1\\rangle \\otimes |0\\rangle \\otimes |1\\rangle \\\\
&= \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\otimes \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} \\otimes \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\\\
&= \\begin{pmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 1 \\end{pmatrix} \\otimes \\begin{pmatrix} 0 \\\\ 1 \\\\ 0 \\\\ 0 \\end{pmatrix} \\\\
&= \\begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\end{pmatrix}^{\\top}

\\end{aligned}
```
Turns out that this effectively result in a lone non zero entry at the index ``13`` when counted from ``0`` to ``2^n``.

```julia-repl
julia> dket(13)
[0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0]

julia> dket(5,n=2)
[0,0,0,0,0,1,0,0]

"""
function dket(x::Int64;n=2)

	if x ==0
		y=[1,0]
	else
		bv=reverse(digits(x,base=2,pad=max(Int(ceil(log2(x))),n)))
    	y=1
		for b in bv
			y = kron(y,ket(b))
		end
	end
	return y
end

"""
Generates the unitary matrix corresponding to a one qubit operation at a specified bit, acting on an ``n-`` register.

A one-qubit gate ``U_2 \\in SU(2)`` acting on ``k``th qubit in an ``n-`` qubit register correspond to a unitary matrix ``U \\in SU(2^n) = \\mathbb{I}_{2}^{\\otimes k} \\otimes U_{2} \\otimes \\mathbb{I}_{2}^{\\otimes n-k}``.

The unitary matrix
```math
U_2 = \\begin{pmatrix} \\acute{a} & \\acute{b} \\\\ -\\acute{b}^{*} & \\acute{a}^{*}\\end{pmatrix} = \\frac{1}{\\lvert a \\rvert^{2}+\\lvert b \\rvert^{2}}\\begin{pmatrix} a & b \\\\ -b^{*} & a^{*}\\end{pmatrix}
```. Here ``a,b \\in \\mathbb{C}, \\lvert \\acute{a} \\rvert^{2}+\\lvert \\acute{b} \\rvert^{2}=1
``
is the embedded single qubit operation.

#### Arguments
* a,b 	- Elements of ``U_2 = \\begin{pmatrix} \\acute{a} & \\acute{b} \\\\ -\\acute{b}^{*} & \\acute{a}^{*}\\end{pmatrix} = \\frac{1}{\\lvert a \\rvert^{2}+\\lvert b \\rvert^{2}}\\begin{pmatrix} a & b \\\\ -b^{*} & a^{*}\\end{pmatrix}``. Here ``a,b \\in \\mathbb{C}, \\lvert \\acute{a} \\rvert^{2}+\\lvert \\acute{b} \\rvert^{2}=1``.
* n     -- Number of registers ``n \\in \\mathbb{N} ``
* k     -- The position of the qubit ``1 \\le k \\le n , k \\in \\mathbb{N} ``
* U    -- Output ``U_2 \\in SU(2^n)``

#### Dependency
Uses the Kronecker.jl package

```julia-repl
julia> oneQubitU(7+3im,1-2im,3,1)
8×8 Matrix{ComplexF64}:
  0.9+0.4im   0.0+0.0im   .0+0.0im   .0+0.0im  …  .0-0.0im  .0-0.0im  0.0-.0im
  0.0+0.0im   0.9+0.4im   .0+0.0im   .0+0.0im     .1-0.3im  .0-0.0im  0.0-.0im
  0.0+0.0im   0.0+0.0im   .9+0.4im   .0+0.0im     .0-0.0im  .1-0.3im  0.0-.0im
  0.0+0.0im   0.0+0.0im   .0+0.0im   .9+0.4im     .0-0.0im  .0-0.0im  0.1-.3im
 -0.1-0.3im  -0.0-0.0im  -.0-0.0im  -.0-0.0im     .0-0.0im  .0-0.0im  0.0-.0im
 -0.0-0.0im  -0.1-0.3im  -.0-0.0im  -.0-0.0im  …  .9-0.4im  .0-0.0im  0.0-.0im
 -0.0-0.0im  -0.0-0.0im  -.1-0.3im  -.0-0.0im     .0-0.0im  .9-0.4im  0.0-.0im
 -0.0-0.0im  -0.0-0.0im  -.0-0.0im  -.1-0.3im     .0-0.0im  .0-0.0im  0.9-.4im

julia> oneQubitU(3,4,3,2)
8×8 Matrix{Float64}:
  0.6   0.0  0.8  0.0   0.0   0.0  0.0  0.0
  0.0   0.6  0.0  0.8   0.0   0.0  0.0  0.0
 -0.8  -0.0  0.6  0.0  -0.0  -0.0  0.0  0.0
 -0.0  -0.8  0.0  0.6  -0.0  -0.0  0.0  0.0
  0.0   0.0  0.0  0.0   0.6   0.0  0.8  0.0
  0.0   0.0  0.0  0.0   0.0   0.6  0.0  0.8
 -0.0  -0.0  0.0  0.0  -0.8  -0.0  0.6  0.0
 -0.0  -0.0  0.0  0.0  -0.0  -0.8  0.0  0.6
```

"""
function oneQubitU(a,b,n,k)
	U0=[a b;-b' a']
	U2= U0./norm([a b])

	L  = (k-1 <= 1) ? I(2) : I(2) ⊗ (k-1)
	R  = (n-k <= 1) ? I(2) : I(2) ⊗ (n-k)
	UL = (k-1 > 0) ? L⊗U2 : U2
	UR = (n-k > 0) ? UL⊗R : UL

	U= Matrix(UR)

end

end