text:
abstract type AbstractBranchResult end

####################################################################################################
namedprintsol(x) = (x = x,)
namedprintsol(x::Real) = (x = x,)
namedprintsol(x::NamedTuple) = x
namedprintsol(x::Tuple) = (;zip((Symbol("x$i") for i in 1:length(x)), x)...)
mergefromuser(x, a::NamedTuple) = merge(namedprintsol(x), a)
####################################################################################################
# Structure to hold continuation result
"""
$(TYPEDEF)

Structure which holds the results after a call to [`continuation`](@ref).

You can see the propertynames of a result by using `propertynames(::ContResult)` or by typing `br.` + TAB where `br::ContResult`.

# Fields

$(TYPEDFIELDS)

# Associated methods
- `length(br)` number of the continuation steps
- `eigenvals(br, ind)` returns the eigenvalues for the ind-th continuation step
- `eigenvec(br, ind, indev)` returns the indev-th eigenvector for the ind-th continuation step
- `br[k+1]` gives information about the k-th step
"""
@with_kw_noshow struct ContResult{Ta, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl <: Lens} <: AbstractBranchResult
	"holds the low-dimensional information about the branch. More precisely, `branch[:, i+1]` contains the following information `(printSolution(u, param), param, itnewton, itlinear, ds, theta, n_unstable, n_imag, stable, step)` for each continuation step `i`.\n
  - `itnewton` number of Newton iterations
  - `itlinear` total number of linear iterations during corrector
  - `n_unstable` number of eigenvalues with positive real part for each continuation step (to detect stationary bifurcation)
  - `n_imag` number of eigenvalues with positive real part and non zero imaginary part for each continuation step (to detect Hopf bifurcation).
  - `stable`  stability of the computed solution for each continuation step. Hence, `stable` should match `eig[step]` which corresponds to `branch[k]` for a given `k`.
  - `step` continuation step (here equal `i`)"
	branch::StructArray{Ta}

	"A vector with eigen-elements at each continuation step."
	eig::Vector{NamedTuple{(:eigenvals, :eigenvec, :step), Tuple{Teigvals, Teigvec, Int64}}}

	"Vector of solutions sampled along the branch. This is set by the argument `saveSolEveryNsteps::Int64` (default 0) in [`ContinuationPar`](@ref)."
	sol::Ts

	"The parameters used for the call to `continuation` which produced this branch."
	contparams::ContinuationPar

	"Type of solutions computed in this branch."
	type::Symbol = :Equilibrium

	"Structure associated to the functional, useful for branch switching. For example, when computing periodic orbits, the functional `PeriodicOrbitTrapProblem`, `ShootingProblem`... will be saved here."
	functional::Tfunc = nothing

	"Parameters passed to continuation and used in the equation `F(x, par) = 0`."
	params::Tpar = nothing

	"Parameter axis used for computing the branch"
	lens::Tl

	"A vector holding the set of detected bifurcation points. See [`SpecialPoint`](@ref) for a description of the fields."
	specialpoint::Vector{Biftype}
end

# returns the number of steps in a branch
Base.length(br::AbstractBranchResult) = length(br.branch)

# check whether the eigenvalues are saved in the branch
# this is a good test bifucause we always fill br.eig with a dummy vector :(
@inline haseigenvalues(br::ContResult) = hasstability(br)
@inline haseigenvalues(br::AbstractBranchResult) = haseigenvalues(br.γ)

# check whether the eigenvectors are saved in the branch
@inline haseigenvector(br::ContResult{Ta, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl} ) where {Ta, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl } = Teigvec != Nothing
@inline haseigenvector(br::AbstractBranchResult) = haseigenvector(br.γ)

@inline hasstability(br::AbstractBranchResult) = typeof(br.branch).parameters[1].parameters[2].parameters[end-1] == Bool

getfirstusertype(br::AbstractBranchResult) = keys(br.branch[1])[1]
@inline getvectortype(br::AbstractBranchResult) = getVectorType(eltype(br.specialpoint))
@inline getvectoreltype(br::AbstractBranchResult) = eltype(getvectortype(br))
setParam(br::AbstractBranchResult, p0) = set(br.params, br.lens, p0)
Base.getindex(br::ContResult, k::Int) = (br.branch[k]..., eigenvals = haseigenvalues(br) ? br.eig[k].eigenvals : nothing, eigenvec = haseigenvector(br) ? br.eig[k].eigenvec : nothing)

function Base.getproperty(br::ContResult, s::Symbol)
	if s in (:specialpoint, :contparams, :lens, :sol, :type, :branch, :eig, :functional, :params)
		getfield(br, s)
	else
		getproperty(br.branch, s)
	end
end
Base.propertynames(br::ContResult) = (propertynames(br.branch)..., :specialpoint, :contparams, :lens, :sol, :type, :branch, :eig, :functional, :params)
@inline kernelDim(br::ContResult, ind) = kernelDim(br.specialpoint[ind])

"""
$(SIGNATURES)

Return the eigenvalues of the ind-th continuation step. `verbose` is used to tell the number of unstable eigen elements.
"""
function eigenvals(br::AbstractBranchResult, ind::Int, verbose::Bool = false)
	@assert br.eig[ind+1].step == ind "Error in indexing eigenvalues. Please open an issue on the website."
	if verbose
		println("--> There are ", br.branch[ind].n_unstable, " unstable eigenvalues")
	println("--> Eigenvalues for continuation step ", br.eig[ind+1].step)
	end
	br.eig[ind+1].eigenvals
end

"""
$(SIGNATURES)

Return the eigenvalues of the ind-th bifurcation point.
"""
eigenvalsfrombif(br::AbstractBranchResult, ind::Int) = br.eig[br.specialpoint[ind].idx].eigenvals

"""
$(SIGNATURES)

Return the indev-th eigenvectors of the ind-th continuation step.
"""
function eigenvec(br::AbstractBranchResult, ind::Int, indev::Int)
	@assert br.eig[ind+1].step == ind "Error in indexing eigenvalues. Please open an issue on the website."
	return geteigenvector(br.contparams.newtonOptions.eigsolver, br.eig[ind+1].eigenvec, indev)
end

function Base.show(io::IO, br::ContResult, comment = "")
	println(io, " ┌─ Branch number of points: ", length(br.branch))
	print(io, " ├─ Branch of ")
	printstyled(io, br.type, comment, color=:light_cyan, bold = true)
	print(io, "\n ├─ Type of vectors: ")
	printstyled(io, getvectortype(br), color=:light_cyan, bold = true)
	print(io, "\n ├─ Parameter ")
	printstyled(io, getLensSymbol(br.lens), color=:light_cyan, bold = true)
	println(io, " from ", br.branch[1].param, " to ", br.branch[end].param,)
	if length(br.specialpoint) > 0
		println(io, " └─ Special points:\n\n (ind_ev = index of the bifurcating eigenvalue e.g. `br.eig[idx].eigenvals[ind_ev]`)\n")
		for ii in eachindex(br.specialpoint)
			_show(io, br.specialpoint[ii], ii, String(getLensSymbol(br.lens)))
		end
	end
end

# this function is important in that it gives the eigenelements corresponding to bp and stored in br. We do not check that bp ∈ br for speed reasons
getEigenelements(br::ContResult{T, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl}, bp::Biftype) where {T, Teigvals, Teigvec, Biftype, Ts, Tfunc, Tpar, Tl} = br.eig[bp.idx]

"""
$(SIGNATURES)
Function is used to initialize the composite type `ContResult` according to the options contained in `contParams`

# Arguments
- `br` result from `getStateSummary`
- `par`: parameters
- `lens`: lens to specify the continuation parameter
- `eiginfo`: eigen-elements (eigvals, eigvecs)
"""
 function _ContResult(printsol, br, x0, par, lens::Lens, eiginfo, contParams::ContinuationPar{T, S, E}, computeEigElements::Bool) where {T, S, E}
	# example of bifurcation point
	bif0 = SpecialPoint(x0, T, namedprintsol(printsol))
	# shall we save full solution?
	sol = contParams.saveSolEveryStep > 0 ? [(x = copy(x0), p = get(par, lens), step = 0)] : nothing
	n_unstable = 0; n_imag = 0; stability = true

	if computeEigElements
		evvectors = contParams.saveEigenvectors ? eiginfo[2] : nothing
		_evvectors = (eigenvals = eiginfo[1], eigenvec = evvectors, step = 0)
	else
		_evvectors = (eigenvals = nothing, eigenvec = nothing, step = 0)
	end
	return ContResult(
		branch = StructArray([br]),
		specialpoint = Vector{typeof(bif0)}(undef, 0),
		eig = computeEigElements ? [_evvectors] : empty([_evvectors]),
		sol = sol,
		contparams =  contParams,
		params = par,
		lens = lens)
end

function _ContResult()
	ContResult(branch = StructArray([(a=1,b=2)]), eig = [(eigenvals = Any[], eigenvec=Any[],step=0)],sol=Any[],contparams = opts_br, lens=(@lens _[1]), bifpoint=Any[]);
end
####################################################################################################
"""
$(TYPEDEF)

A Branch is a structure which encapsulates the result of the computation of a branch bifurcating from a bifurcation point.

$(TYPEDFIELDS)
"""
struct Branch{T <: Union{ContResult, Vector{ContResult}}, Tbp} <: AbstractBranchResult
	"Set of branches branching  off the bifurcation point `bp`"
	γ::T
	"Bifurcation point. It is thought as the root of the branches in γ"
	bp::Tbp
end

Base.length(br::Branch) = length(br.γ)

"""
$(SIGNATURES)

Return the bifurcation point of a `::Branch`.
"""
from(br::Branch) = br.bp
from(br::Vector{Branch}) = length(br) > 0 ? from(br[1]) : nothing
from(tree::ContResult) = nothing
getfirstusertype(br::Branch) = getfirstusertype(br.γ)
Base.show(io::IO, br::Branch{T, Tbp}) where {T <: ContResult, Tbp} = show(io, br.γ, " from $(type(br.bp)) bifurcation point.")

# extend the getproperty for easy manipulation of a Branch
# for example, it allows to use the plot recipe for ContResult as is
Base.getproperty(br::Branch, s::Symbol) = s in (:γ, :bp) ? getfield(br, s) : getproperty(br.γ, s)
Base.propertynames(br::Branch) = ((:γ, :bp)..., propertynames(br.γ)...)
Base.getindex(br::Branch, k::Int) = getindex(br.γ, k)
####################################################################################################
_reverse!(x) = reverse!(x)
_reverse!(::Nothing) = nothing
function _reverse(br0::ContResult)
	br = deepcopy(br0)
	nb = length(br.branch)
	if ~isnothing(br.branch)
		@set! br.branch =
			StructArray([setproperties(pt; step = nb - pt.step - 1) for pt in Iterators.reverse(br.branch)])
	end

	if ~isnothing(br.specialpoint)
		@set! br.specialpoint =
			[setproperties(pt;
				step = nb - pt.step - 1,
				idx = nb - pt.idx + 1,
				δ = (-pt.δ[1], -pt.δ[2])) for pt in Iterators.reverse(br.specialpoint)]
	end

	if ~isnothing(br.eig)
		@set! br.eig =
			[setproperties(pt; step = nb - pt.step - 1) for pt in Iterators.reverse(br.eig)]
	end

	if ~isnothing(br.sol)
		@set! br.sol =
			[setproperties(pt; step = nb - pt.step - 1) for pt in Iterators.reverse(br.sol)]
	end
	return br
end

_append!(x,y) = append!(x,y)
_append!(x,::Nothing) = nothing

"""
$(SIGNATURES)

Merge two `ContResult`s and put the result in `br`.
"""
function _cat!(br::ContResult, br2::ContResult)
	# br = deepcopy(br1)
	nb = length(br.branch)
	if ~isnothing(br.branch)
		append!(br.branch,
			[setproperties(pt; step = nb + pt.step) for pt in br2.branch])
	end
	if ~isnothing(br.specialpoint)
		append!(br.specialpoint,
			[setproperties(pt;
				step = nb + pt.step,
				idx = nb + pt.idx) for pt in br2.specialpoint])
	end

	if ~isnothing(br.eig)
		append!(br.eig,
			[setproperties(pt; step = nb + pt.step) for pt in br2.eig])
	end

	if ~isnothing(br.sol)
		append!(br.sol,
			[setproperties(pt; step = nb + pt.step) for pt in br2.sol])
	end
	return br
end

# _catrev(br1::ContResult, br2::ContResult) = _merge!(_reverse(br1), br2)
# _cat(br1::ContResult, br2::ContResult) = _merge!(deepcopy(br1), br2)

"""
Same as _cat! but determine the ordering so that the branches merge properly
"""
function _merge(br1::ContResult, br2::ContResult; tol = 1e-6)
	# find the intersection point
	dst(x1,p1,x2,p2) = max(abs(x1-x2),abs(p1-p2))
	dst(i,j) = dst(br1.branch[i][1],br1.branch[i].param,br2.branch[j][1],br2.branch[j].param)
	ind = (1,1)
	for i in [1,length(br1)], j in [1,length(br2)]
		if dst(i,j) < tol
			ind = (i,j)
			break
		end
	end

	if ind[1] == 1
		if ind[2] == 1
			return _cat!(_reverse(br2),br1)
		else
			return _cat!((br2),br1)
		end
	else
		if ind[2] == 1
			return _cat!(br1,br2)
		else
			return _cat!(br1,_reverse(br2))
		end

	end

	if minimum(br1.branch.param) < minimum(br2.branch.param)
		@debug "b1-b2"
		# br1 is the first branch and then br2
		# we need to look at the indexing
		return _cat!(_reverse(br1), br2)
	else
		@debug "b2-b1"
		return _cat!(_reverse(br1), br2)
	end
end


text:
"""
Main module for `LLLplus.jl` -- lattice reduction and related tools for Julia.

As an example of the functions in the package, see [`lll`](@ref), which does
Lenstra–Lenstra–Lovász lattice reduction of a matrix.

"""
module LLLplus

using LinearAlgebra
using Printf

export
    lll,
    cvp,
    svp,
    brun,
    gauss,
    sizereduction,
    seysen,
    vblast,
    subsetsum,mdsubsetsum,
    integerfeasibility,
    rationalapprox,
    hardsphere, hard_sphere,
    issizereduced,islllreduced,orthogonalitydefect,hermitefactor,seysencond

include("lll.jl")          # lll, gauss, sizereduction
include("cvp.jl")          # cvp, svp
include("brun.jl")
include("seysen.jl")
include("vblast.jl")
include("applications.jl") # subsetsum, integerfeasibility, rationalapprox
include("utilities.jl")

include("hard_sphere.jl")  # may be deprecated in future

roundf(r::Td) where {Td<:Complex} = round(real(r)) + im*round(imag(r));
roundf(r) = round(r);

end # LLLplus module

ext: jl


text:
using NaturalSort
using Test

# write your own tests here
STRS = [
    "",
    "a#",
    "a00",
    "a0a",
    "a1",
    "a1b1",
    "a1b2",
    "a01b3",
    "a1b11",
    "a2",
    "a11",
    "a11b1",
    "a11b2",
    "a11b11",
    "a12",
    "a12b1",
    "a22",
    "abc",
    "b"
]

for i = 1:length(STRS)
    for j = 1:i-1
        @test natural(STRS[j], STRS[i])
    end
    for j = i:length(STRS)
        @test !natural(STRS[j], STRS[i])
    end
end

EQ_STRS = [
    ("a01", "a1"),
    ("a01b2", "a1b2"),
    ("a0", "a00")
]
for (str, eqto) in EQ_STRS
    for i = 1:findfirst(x->x == eqto, STRS)-1
        @test natural(STRS[i], str)
        @test !natural(str, STRS[i])
    end
    for i = findfirst(x->x == eqto, STRS)+1:length(STRS)
        @test !natural(STRS[i], str)
        @test natural(str, STRS[i])
    end
end


text:
import ZXCalculus: ZXDiagram

function ZXDiagram(c::AbstractBlock)
	zxd = ZXDiagram(nqubits(c))
	push_gate!(zxd, c)
	return zxd
end

function push_gate!(zxd::ZXDiagram, c::AbstractBlock)
	push_gate!(zxd, decompose_zx(c))
end

# rotation blocks
function push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,RotationGate{1,T,XGate}}) where {N,T}
	push_gate!(zxd, Val(:X), c.locs[1], c.content.theta/π)
end
function push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,RotationGate{1,T,ZGate}}) where {N,T}
	push_gate!(zxd, Val(:Z), c.locs[1], c.content.theta/π)
end

function push_gate!(zxd::ZXDiagram, c::ChainBlock{N}) where {N}
	push_gate!.(Ref(zxd), subblocks(c))
	zxd
end

# constant block
function push_gate!(zxd::ZXDiagram, c::PutBlock{N,1,HGate}) where {N}
	push_gate!(zxd, Val(:H), c.locs[1])
end

# control blocks
function push_gate!(zxd::ZXDiagram, c::ControlBlock{N,XGate,1}) where {N}
	cloc = c.ctrl_locs[1]
	if c.ctrl_config[1] == 1
		push_gate!(zxd, Val(:CNOT), c.locs[1], cloc)
	else
		push_gate!(zxd, Val(:X), cloc, 1//1)
		push_gate!(zxd, Val(:CNOT), c.locs[1], cloc)
		push_gate!(zxd, Val(:X), cloc, 1//1)
	end
end
function push_gate!(zxd::ZXDiagram, c::ControlBlock{N,ZGate,1}) where {N}
	cloc = c.ctrl_locs[1]
	if c.ctrl_config[1] == 1
		push_gate!(zxd, Val(:CZ), cloc, c.locs[1])
	else
		push_gate!(zxd, Val(:X), cloc, 1//1)
		push_gate!(zxd, Val(:CZ), c.locs[1], cloc)
		push_gate!(zxd, Val(:X), cloc, 1//1)
	end
end

function push_gate!(zxd::ZXDiagram, c::PutBlock{N,2,SWAPGate}) where {N}
	a, b = c.locs
	push_gate!(zxd, Val(:SWAP), [a, b])
end

# ref: https://qiskit.org/textbook/ch-gates/more-circuit-identities.html
function decompose_zx(c::ControlBlock{N,XGate,2}) where N
	a, b = getclocs(c)
	loc = c.locs[1]
	chain(N,
		put(loc=>H),
		cnot(b, loc),
		put(loc=>ConstGate.Tdag),
		cnot(a, loc),
		put(loc=>ConstGate.T),
		cnot(b, loc),
		put(loc=>ConstGate.Tdag),
		cnot(a, loc),
		put(loc=>ConstGate.T),
		put(b=>ConstGate.T),
		put(loc=>H),
		cnot(a, b),
		put(a=>ConstGate.T),
		put(b=>ConstGate.Tdag),
		cnot(a, b),
	)
end

function decompose_zx(c::ControlBlock{N,YGate,1}) where {N}
	a = getclocs(c)[1]
	loc = c.locs[1]
	chain(N, put(loc=>ConstGate.Sdag), cnot(a, loc), put(loc=>ConstGate.S))
end

function decompose_zx(c::ControlBlock{N,ShiftGate{T},1}) where {N,T}
	a = getclocs(c)[1]
	loc = c.locs[1]
	θ = c.content.theta
	chain(N, put(a=>ShiftGate(θ/2)), put(loc=>Rz(θ/2)), cnot(a,loc), put(loc=>Rz(-θ/2)), cnot(a,loc))
end

function decompose_zx(c::ControlBlock{N,RotationGate{1,T,YGate},1}) where {N,T}
	a = getclocs(c)[1]
	loc = c.locs[1]
	θ = c.content.theta
	chain(N, put(loc=>Ry(θ/2)), cnot(a,loc), put(loc=>Ry(-θ/2)), cnot(a,loc))
end

function decompose_zx(c::ControlBlock{N,RotationGate{1,T,XGate},1}) where {N,T}
	a = getclocs(c)[1]
	loc = c.locs[1]
	θ = c.content.theta
	chain(N, put(loc=>Rx(θ/2)), cz(a,loc), put(loc=>Rx(-θ/2)), cz(a,loc))
end

function decompose_zx(c::ControlBlock{N,RotationGate{1,T,ZGate},1}) where {N,T}
	a = getclocs(c)[1]
	loc = c.locs[1]
	θ = c.content.theta
	chain(N, put(loc=>Rz(θ/2)), cnot(a,loc), put(loc=>Rz(-θ/2)), cnot(a,loc))
end

# constant block
function decompose_zx(c::PutBlock{N,1,RotationGate{1,T,XGate}}) where {N,T}
	put(N, c.locs[1]=>Rx(π))
end
function decompose_zx(c::PutBlock{N,1,RotationGate{1,T,ZGate}}) where {N,T}
	put(N, c.locs[1]=>Rz(π))
end
function decompose_zx(c::PutBlock{N,1,ConstGate.SGate}) where {N}
	put(N, c.locs[1]=>Rz(π/2))
end
function decompose_zx(c::PutBlock{N,1,ConstGate.TGate}) where {N}
	put(N, c.locs[1]=>Rz(π/4))
end
function decompose_zx(c::PutBlock{N,1,ConstGate.SdagGate}) where {N}
	put(N, c.locs[1]=>Rz(-π/2))
end
function decompose_zx(c::PutBlock{N,1,ConstGate.TdagGate}) where {N}
	put(N, c.locs[1]=>Rz(-π/4))
end
function decompose_zx(c::PutBlock{N,1,ShiftGate{T}}) where {N,T}
	put(N, c.locs[1]=>Rz(c.content.theta))
end

function getclocs(c::ControlBlock)
	(2 .* c.ctrl_config .- 1) .* c.ctrl_locs
end


text:
using Makie, GLFW, GeometryTypes, Reactive, FileIO, ColorBrewer, Colors
using GLVisualize
using GLVisualize: loadasset, assetpath

function xy_data(x, y)
    r = sqrt(x*x + y*y)
    r == 0.0 ? 1f0 : (sin(r)/r)
end

function custom_theme(scene)
    @theme theme = begin
        linewidth = to_float(3)
        colormap = to_colormap(:RdYlGn)#to_colormap(:RdPu)
        scatter = begin
            marker = to_spritemarker(Circle)
            markersize = to_float(0.03)
            strokecolor = to_color(:white)
            strokewidth = to_float(0.01)
            glowcolor = to_color(RGBA(0, 0, 0, 0.4))
            glowwidth = to_float(0.1)
        end
    end
    # update theme values
    scene[:theme] = theme
end

function testfunc_1()

    img = loadasset("doge.png")
    scene = Scene(resolution = (500, 500))

    display(scene)
    show(scene)
    println(scene)

    is = image(img)
    center!(scene)
    subscene = Scene(scene, Signal(SimpleRectangle(0, 0, 200, 200)))
    scatter(subscene, rand(100) * 200, rand(100) * 200, markersize = 4)
    center!(scene)

end
function testfunc_2()

    scene = Scene(resolution = (500, 500));
    x = [0, 1, 2, 0];
    y = [0, 0, 1, 2];
    z = [0, 2, 0, 1];
    color = [:red, :green, :blue, :yellow];
    i = [0, 0, 0, 1];
    j = [1, 2, 3, 2];
    k = [2, 3, 1, 3];

    indices = [1, 2, 3, 1, 3, 4, 1, 4, 2, 2, 3, 4];
    mesh(x, y, z, indices, color = color);
    r = linspace(-0.5, 2.5, 4);
    axis(r, r, r);
    center!(scene);

end
function testfunc_3()

    scene = Scene(resolution = (500, 500))
    Makie.Makie.volume(rand(32, 32, 32), algorithm = :iso)
    center!(scene)

end
function testfunc_4()

    scene = Scene(resolution = (500, 500))
    heatmap(rand(32, 32))
    center!(scene)

end
function testfunc_5()

    scene = Scene(resolution = (500, 500))
    r = linspace(-10, 10, 512)
    z = ((x, y)-> sin(x) + cos(y)).(r, r')
    Makie.contour(r, r, z, levels = 5, color = ColorBrewer.palette("RdYlBu", 5))
    center!(scene)

end
function testfunc_6()

    scene = Scene(resolution = (500, 500))
    vx = -1:0.1:1;
    vy = -1:0.1:1;

    f(x, y) = (sin(x*10) + cos(y*10)) / 4
    psurf = surface(vx, vy, f)

    pos = lift_node(psurf[:x], psurf[:y], psurf[:z]) do x, y, z
        vec(Point3f0.(x, y', z .+ 0.5))
    end
    pscat = scatter(pos)
    plines = lines(view(pos, 1:2:length(pos)))
    center!(scene)
    @theme theme = begin
        markersize = to_markersize2d(0.01)
        strokecolor = to_color(:white)
        strokewidth = to_float(0.01)
    end
    # this pushes all the values from theme to the plot
    push!(pscat, theme)
    pscat[:glow_color] = to_node(RGBA(0, 0, 0, 0.4), x->to_color((), x))
    # apply it to the scene
    custom_theme(scene)
    # From now everything will be plotted with new theme
    psurf = surface(vx, 1:0.1:2, psurf[:z])
    center!(scene)

end
function testfunc_7()

    scene = Scene(resolution = (500, 500))
    sv = scatter(rand(Point3f0, 100))
    similar(sv, rand(10), rand(10), rand(10), color = :black, markersize = 0.4)
    scene

end
function testfunc_8()

    scene = Scene(resolution = (500, 500))
    large_sphere = HyperSphere(Point3f0(0), 1f0)
    positions = decompose(Point3f0, large_sphere)
    colS = [Colors.RGBA{Float32}(rand(), rand(), rand(), 1.) for i = 1:length(positions)]
    sizesS = [rand(Vec3f0) .* 0.5f0 for i = 1:length(positions)]
    meshscatter(positions, color = colS, markersize = sizesS)
    center!(scene)

end
function testfunc_9()

    scene = Scene(resolution = (500, 500))
    y = [
        -0.997669
        -0.979084
        -0.942261
        -0.887885
        -0.81697
        -0.730836
        -0.631088
        -0.519584
        -0.398401
        -0.269797
        -0.136167
        0.0
        0.136167
        0.269797
        0.398401
        0.519584
        0.631088
        0.730836
        0.81697
        0.887885
        0.942261
        0.979084
        0.997669
    ]
    contour(linspace(-0.99, 0.99, 23), y, rand(23, 23), levels = 10)
    center!(scene)

end
function testfunc_10()

    scene = Scene(resolution = (500, 500))
    # define points/edges
    perturbfactor = 4e1
    N = 3; nbfacese = 30; radius = 0.02
    large_sphere = HyperSphere(Point3f0(0), 1f0)
    positions = decompose(Point3f0, large_sphere, 30)
    np = length(positions)
    pts = [positions[k][l] for k = 1:length(positions), l = 1:3]
    pts = vcat(pts, 1.1 * pts + randn(size(pts)) / perturbfactor) # light position influence ?
    edges = hcat(collect(1:np), collect(1:np) + np)
    ne = size(edges, 1); np = size(pts, 1)
    # define markers meshes
    meshC = GeometryTypes.GLNormalMesh(GeometryTypes.Cylinder{3, Float32}(
                                       GeometryTypes.Point3f0(0., 0., 0.),
                                       GeometryTypes.Point3f0(0., 0, 1.),
                                       Float32(1)), nbfacese)

    meshS = GeometryTypes.GLNormalMesh(large_sphere, 20)
    # define colors, markersizes and rotations
    pG = [GeometryTypes.Point3f0(pts[k, 1], pts[k, 2], pts[k, 3]) for k = 1:np]
    lengthsC = sqrt.(sum((pts[edges[:,1], :] .- pts[edges[:, 2], :]) .^ 2, 2))
    sizesC = [GeometryTypes.Vec3f0(radius, radius, lengthsC[i]) for i = 1:ne]
    sizesC = [Vec3f0(1., 1., 1.) for i = 1:ne]
    colorsp = [Colors.RGBA{Float32}(rand(), rand(), rand(), 1.) for i = 1:np]
    colorsC = [(colorsp[edges[i, 1]] + colorsp[edges[i, 2]]) / 2. for i = 1:ne]
    sizesC = [Vec3f0(radius, radius, lengthsC[i]) for i = 1:ne]
    Qlist = zeros(ne, 4)
    for k = 1:ne
        ct = GeometryTypes.Cylinder{3, Float32}(
                    GeometryTypes.Point3f0(pts[edges[k, 1], 1], pts[edges[k, 1], 2], pts[edges[k, 1], 3]),
                    GeometryTypes.Point3f0(pts[edges[k, 2], 1], pts[edges[k, 2], 2], pts[edges[k, 2], 3]),
                    Float32(1))
        Q = GeometryTypes.rotation(ct)
        r = 0.5 * sqrt(1 + Q[1, 1] + Q[2, 2] + Q[3, 3]); Qlist[k, 4] = r
        Qlist[k, 1] = (Q[3, 2] - Q[2, 3]) / (4 * r)
        Qlist[k, 2] = (Q[1, 3] - Q[3, 1]) / (4 * r)
        Qlist[k, 3] = (Q[2, 1] - Q[1, 2]) / (4 * r)
    end
    rotationsC = AbstractVector[Vec4f0(Qlist[i, 1], Qlist[i, 2], Qlist[i, 3], Qlist[i, 4]) for i = 1:ne]
    # plot
    hm = Makie.meshscatter(pG[edges[:, 1]], color = colorsC, marker = meshC,
                           markersize = sizesC,  rotations = rotationsC)
    hp = Makie.meshscatter(pG, color = colorsp, marker = meshS, markersize = radius)

    r = linspace(-1.3, 1.3, 4); Makie.axis(r, r, r)
    center!(scene)

end
function testfunc_11()

    scene = Scene(resolution = (500, 500))
    large_sphere = HyperSphere(Point3f0(0), 1f0)
    positions = decompose(Point3f0, large_sphere)
    linepos = view(positions, rand(1:length(positions), 1000))
    lines(linepos, linewidth = 0.1, color = :black)
    scatter(positions, strokewidth = 0.02, strokecolor = :white, color = RGBA(0.9, 0.2, 0.4, 0.6))
    r = linspace(-1.5, 1.5, 5)
    axis(r, r, r)
    scene

end
function testfunc_12()

    scene = Scene(resolution = (500, 500))
    large_sphere = HyperSphere(Point3f0(0), 1f0)
    positions = decompose(Point3f0, large_sphere)
    meshscatter(positions, color = RGBA(0.9, 0.2, 0.4, 1))
    scene

end
function testfunc_13()

    scene = Scene(resolution = (500, 500))

    r = linspace(-2, 2, 40)
    surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
    z = surf_func(20)
    surf = surface(r, r, z)

    wf = wireframe(r, r, surf[:z] .+ 1.0,
        linewidth = 2f0, color = lift_node(x-> x[5], surf[:colormap])
    )
    xy = linspace(-2.1, 2.1, 4)
    axis(xy, xy, linspace(0, 2, 4))
    center!(scene)

    io = VideoStream(scene)
    for i in linspace(0, 60, 100)
        surf[:z] = surf_func(i)
        recordframe!(io)
    end
    scene

end
function testfunc_14()

    scene = Scene(resolution = (500, 500))

    N = 40
    r = linspace(-2, 2, 40)
    surf_func(i) = [Float32(xy_data(x*i, y*i)) for x = r, y = r]
    surface(
        r, r, surf_func(10),
        color = GLVisualize.loadasset("doge.png")
    )
    center!(scene)
    scene

end
function testfunc_15()

    scene = Scene(resolution = (500, 500))
    x = GLVisualize.loadasset("cat.obj")
    Makie.mesh(x.vertices, x.faces, color = :black)
    pos = map(x.vertices, x.normals) do p, n
        p => p .+ (normalize(n) .* 0.05f0)
    end
    linesegment(pos)
    scene

end
function testfunc_16()

    scene = Scene(resolution = (500, 500))
    mesh(GLVisualize.loadasset("cat.obj"))
    r = linspace(-0.1, 1, 4)
    center!(scene)
    scene

end
function testfunc_17()

    scene = Scene(resolution = (500, 500))
    cat = load(assetpath("cat.obj"), GLNormalUVMesh)
    Makie.mesh(cat, color = loadasset("diffusemap.tga"))
    center!(scene)

end
function testfunc_18()

    scene = Scene(resolution = (500, 500))
    Makie.mesh(Sphere(Point3f0(0), 1f0))
    center!(scene)
    scene

end
function testfunc_19()

    scene = Scene(resolution = (500, 500))
    wireframe(GLVisualize.loadasset("cat.obj"))
    center!(scene)
    scene

end
function testfunc_20()

    scene = Scene(resolution = (500, 500))
    wireframe(Sphere(Point3f0(0), 1f0))
    center!(scene)
    scene

end
function testfunc_21()

    scene = Scene(resolution = (500, 500))
    heatmap(rand(32, 32))
    center!(scene)

end
function testfunc_22()

    scene = Scene(resolution = (500, 500), color = :black)
    stars = 100_000
    scatter((rand(Point3f0, stars) .- 0.5) .* 10,
        glowwidth = 0.005, glow_color = :white, color = RGBA(0.8, 0.9, 0.95, 0.4),
        markersize = rand(linspace(0.0001, 0.01, 100), stars)
    )
    scene

end
function testfunc_23()

    scene = Scene(resolution = (500, 500))
    Makie.volume(rand(32, 32, 32), algorithm = :iso)
    center!(scene)

end
function testfunc_24()

    scene = Scene(resolution = (500, 500))
    scatter(Point3f0[(1,0,0), (0,1,0), (0,0,1)], marker=[:x, :circle, :cross])
    axis(scene, linspace(0, 1, 4), linspace(0, 1, 4), linspace(0, 1, 4))
    center!(scene);

end
function testfunc_25()

    scene = Scene(resolution = (500, 500))

    x = map([:dot, :dash, :dashdot], [2, 3, 4]) do ls, lw
        linesegment(linspace(1, 5, 100), rand(100), rand(100), linestyle = ls, linewidth = lw)
    end
    push!(x, scatter(linspace(1, 5, 100), rand(100), rand(100)))
    center!(scene)
    l = Makie.legend(x, ["attribute $i" for i in 1:4])
    l[:position] = (0, 1)
    l[:backgroundcolor] = RGBA(0.95, 0.95, 0.95)
    l[:strokecolor] = RGB(0.8, 0.8, 0.8)
    l[:gap] = 30
    l[:textsize] = 19
    l[:linepattern] = Point2f0[(0,-0.2), (0.5, 0.2), (0.5, 0.2), (1.0, -0.2)]
    l[:scatterpattern] = decompose(Point2f0, Circle(Point2f0(0.5, 0), 0.3f0), 9)
    l[:markersize] = 2f0
    scene

end
function testfunc_26()

    scene = Scene(resolution = (500, 500))
    cmap = collect(linspace(to_color(:red), to_color(:blue), 20))
    l = Makie.legend(cmap, 1:4)
    l[:position] = (1.0,1.0)
    l[:textcolor] = :blue
    l[:strokecolor] = :black
    l[:strokewidth] = 1
    l[:textsize] = 15
    l[:textgap] = 5
    scene
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

include("ZMQCM.jl")

# @spawn run(`julia broker.jl`)

start_master(parse(Int,ARGS[1]))

resp = pmap(x -> myid() *2, [1:nworkers()])

println(resp)

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
using GeometryBasics
using Statistics

@cell "Test heatmap + image overlap" begin
    heatmap(RNG.rand(32, 32))
    image!(map(x -> RGBAf0(x, 0.5, 0.5, 0.8), RNG.rand(32, 32)))
end

@cell "poly and colormap" begin
    # example by @Paulms from JuliaPlots/Makie.jl#310
    points = Point2f0[[0.0, 0.0], [0.1, 0.0], [0.1, 0.1], [0.0, 0.1]]
    colors = [0.0 ,0.0, 0.5, 0.0]
    scene = poly(points, color=colors, colorrange=(0.0, 1.0))
    points = Point2f0[[0.1, 0.1], [0.2, 0.1], [0.2, 0.2], [0.1, 0.2]]
    colors = [0.5,0.5,1.0,0.3]
    poly!(scene, points, color=colors, colorrange=(0.0, 1.0))
    scene
end

@cell "quiver" begin
    x = range(-2, stop=2, length=21)
    arrows(x, x, RNG.rand(21, 21), RNG.rand(21, 21), arrowsize=0.05)
end

@cell "Arrows on hemisphere" begin
    s = Sphere(Point3f0(0), 0.9f0)
    scene = mesh(s, transparency=true, alpha=0.05)
    pos = decompose(Point3f0, s)
    dirs = decompose_normals(s)
    arrows!(scene, pos, dirs, arrowcolor=:red, arrowsize=0.1, linecolor=:red)
end

@cell "image" begin
    vbox(
        image(AbstractPlotting.logo(), scale_plot=false),
        image(RNG.rand(100, 500), scale_plot=false),
    )
end

@cell "FEM polygon 2D" begin
    coordinates = [
        0.0 0.0;
        0.5 0.0;
        1.0 0.0;
        0.0 0.5;
        0.5 0.5;
        1.0 0.5;
        0.0 1.0;
        0.5 1.0;
        1.0 1.0;
    ]
    connectivity = [
        1 2 5;
        1 4 5;
        2 3 6;
        2 5 6;
        4 5 8;
        4 7 8;
        5 6 9;
        5 8 9;
    ]
    color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
    poly(coordinates, connectivity, color=color, strokecolor=(:black, 0.6), strokewidth=4)
end

@cell "FEM mesh 2D" begin
    coordinates = [
        0.0 0.0;
        0.5 0.0;
        1.0 0.0;
        0.0 0.5;
        0.5 0.5;
        1.0 0.5;
        0.0 1.0;
        0.5 1.0;
        1.0 1.0;
    ]
    connectivity = [
        1 2 5;
        1 4 5;
        2 3 6;
        2 5 6;
        4 5 8;
        4 7 8;
        5 6 9;
        5 8 9;
    ]
    color = [0.0, 0.0, 0.0, 0.0, -0.375, 0.0, 0.0, 0.0, 0.0]
    scene = mesh(coordinates, connectivity, color=color, shading=false)
    wireframe!(scene[end][1], color=(:black, 0.6), linewidth=3)
end

@cell "colored triangle" begin
    mesh(
        [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)], color=[:red, :green, :blue],
        shading=false
    )
end

@cell "colored triangle" begin
    poly(
        [(0.0, 0.0), (0.5, 1.0), (1.0, 0.0)],
        color=[:red, :green, :blue],
        strokecolor=:black, strokewidth=2
    )
end

@cell "Subscenes" begin
    img = RNG.rand(RGBAf0, 100, 100)
    scene = image(img, show_axis=false)
    subscene = Scene(scene, IRect(100, 100, 300, 300))
    scatter!(subscene, RNG.rand(100) * 200, RNG.rand(100) * 200, markersize=4)
    scene
end

@cell "scale_plot" begin
    t = range(0, stop=1, length=500) # time steps
    θ = (6π) .* t    # angles
    x =  # x coords of spiral
    y =  # y coords of spiral
    p1 = lines(t .* cos.(θ), t .* sin.(θ);
        color=t, colormap=:algae, linewidth=8,
        scale_plot=false)
end

@cell "Polygons" begin
    scene = Scene(resolution=(500, 500))
    points = decompose(Point2f0, Circle(Point2f0(50), 50f0))
    pol = poly!(scene, points, color=:gray, strokewidth=10, strokecolor=:red)
    # Optimized forms
    poly!(scene, [Circle(Point2f0(50 + 300), 50f0)], color=:gray, strokewidth=10, strokecolor=:red)
    poly!(scene, [Circle(Point2f0(50 + i, 50 + i), 10f0) for i = 1:100:400], color=:red)
    poly!(scene, [FRect2D(50 + i, 50 + i, 20, 20) for i = 1:100:400], strokewidth=2, strokecolor=:green)
    linesegments!(scene,
        [Point2f0(50 + i, 50 + i) => Point2f0(i + 70, i + 70) for i = 1:100:400], linewidth=8, color=:purple
    )
end

@cell "Contour Function" begin
    r = range(-10, stop=10, length=512)
    z = ((x, y) -> sin(x) + cos(y)).(r, r')
    contour(r, r, z, levels=5, colormap=:viridis, linewidth=3)
end

@cell "Hbox" begin
    t = range(-122277.9, stop=-14798.0, length=29542)
    x = -42 .- RNG.randn(length(t))
    sc1 = scatter(t, x, color=:black, markersize=1.0)
    sc2 = lines(t[1:end - 1], diff(x), color=:blue)
    hbox(sc2, sc1)
end

@cell "Customize Axes" begin
    x = LinRange(0, 3pi, 200); y = sin.(x)
    lin = lines(x, y, padding=(0.0, 0.0), axis=(names = (axisnames = ("", ""),),
        grid = (linewidth = (0, 0),),))
end

@cell "contour" begin
    y = range(-0.997669, stop=0.997669, length=23)
    contour(range(-0.99, stop=0.99, length=23), y, RNG.rand(23, 23), levels=10)
end

@cell "Text Annotation" begin
    text(
        ". This is an annotation!",
        position=(300, 200),
        align=(:center,  :center),
        textsize=60,
        font="Blackchancery"
    )
end

@cell "Text rotation" begin
    scene = Scene()
    pos = (500, 500)
    posis = Point2f0[]
    for r in range(0, stop=2pi, length=20)
        p = pos .+ (sin(r) * 100.0, cos(r) * 100)
        push!(posis, p)
        t = text!(
            scene, "test",
            position=p,
            textsize=50,
            rotation=1.5pi - r,
            align=(:center, :center)
        )
    end
    scatter!(scene, posis, markersize=10)
end

@cell "Standard deviation band" begin
    # Sample 100 Brownian motion path and plot the mean trajectory together
    # with a ±1σ band (visualizing uncertainty as marginal standard deviation).
    n, m = 100, 101
    t = range(0, 1, length=m)
    X = cumsum(RNG.randn(n, m), dims=2)
    X = X .- X[:, 1]
    μ = vec(mean(X, dims=1)) # mean
    lines(t, μ)              # plot mean line
    σ = vec(std(X, dims=1))  # stddev
    band!(t, μ + σ, μ - σ)   # plot stddev band
end

@cell "Streamplot animation" begin
    v(x::Point2{T}, t) where T = Point2{T}(one(T) * x[2] * t, 4 * x[1])
    sf = Node(Base.Fix2(v, 0e0))
    title_str = Node("t = 0.00")
    sp = streamplot(sf, -2..2, -2..2;
                    linewidth=2, padding=(0, 0),
                    arrow_size=0.09, colormap=:magma)
    sc = title(sp, title_str)
    Record(sc, LinRange(0, 20, 5)) do i
        sf[] = Base.Fix2(v, i)
        title_str[] = "t = $(round(i; sigdigits=2))"
    end
end

@cell "Line changing colour" begin
    scene = lines(RNG.rand(10); linewidth=10)

    Record(scene, 1:255; framerate=60) do i
        scene.plots[2][:color] = RGBf0(i / 255, (255 - i) / 255, 0) # animate scene
    end
end

let
    struct FitzhughNagumo2
    end
    (()-> FitzhughNagumo2())()
end

let
    struct FitzhughNagumo{T}
        ϵ::T
        s::T
        γ::T
        β::T
    end
    @cell "streamplot" begin
        P = FitzhughNagumo(0.1, 0.0, 1.5, 0.8)
        f(x, P::FitzhughNagumo) = Point2f0(
            (x[1] - x[2] - x[1]^3 + P.s) / P.ϵ,
            P.γ * x[1] - x[2] + P.β
        )
        f(x) = f(x, P)
        streamplot(f, -1.5..1.5, -1.5..1.5, colormap=:magma)
    end
end

@cell "Transforming lines" begin
    N = 7 # number of colours in default palette
    sc = Scene()
    st = Stepper(sc)

    xs = 0:9        # data
    ys = zeros(10)

    for i in 1:N    # plot lines
        lines!(sc,
            xs, ys;
            color=AbstractPlotting.default_palettes.color[][i],
            limits=FRect((0, 0), (10, 10)),
            linewidth=5
        ) # plot lines with colors
    end

    center!(sc)

    step!(st)

    for (i, rot) in enumerate(LinRange(0, π / 2, N))
        AbstractPlotting.rotate!(sc.plots[i + 1], rot)
        arc!(sc,
            Point2f0(0),
            (8 - i),
            pi / 2,
            (pi / 2 - rot);
            color=sc.plots[i + 1].color,
            linewidth=5,
            linestyle=:dash
        )
    end

    step!(st)
end

@cell "Errorbars x y low high" begin
    x = 1:10
    y = sin.(x)
    scene = scatter(x, y)
    errorbars!(scene, x, y, RNG.rand(10) .+ 0.5, RNG.rand(10) .+ 0.5)
    errorbars!(scene, x, y, RNG.rand(10) .* 0.3 .+ 0.1, RNG.rand(10) .* 0.3 .+ 0.1,
        color=:red, direction=:x)
    scene
end

@cell "Errorbars xy low high" begin
    x = 1:10
    y = sin.(x)
    xy = Point2f0.(x, y)
    scene = scatter(xy)
    errorbars!(scene, xy, RNG.rand(10) .+ 0.5, RNG.rand(10) .+ 0.5)
    errorbars!(scene, xy, RNG.rand(10) .* 0.3 .+ 0.1, RNG.rand(10) .* 0.3 .+ 0.1,
        color=:red, direction=:x)
    scene
end

@cell "Errorbars xy error" begin
    x = 1:10
    y = sin.(x)
    xy = Point2f0.(x, y)
    scene = scatter(xy)
    errorbars!(scene, xy, RNG.rand(10) .+ 0.5)
    errorbars!(scene, xy, RNG.rand(10) .* 0.3 .+ 0.1,
        color=:red, direction=:x)
    scene
end

@cell "Simple pie chart" begin
    scene, layout = layoutscene(resolution=(800, 800))
    ax = layout[1, 1] = LAxis(scene, autolimitaspect=1)

    pie!(ax, 1:5, color=1:5)

    scene
end

@cell "Hollow pie chart" begin
    scene, layout = layoutscene(resolution=(800, 800))
    ax = layout[1, 1] = LAxis(scene, autolimitaspect=1)

    pie!(ax, 1:5, color=1:5, radius=2, inner_radius=1)

    scene
end

@cell "Open pie chart" begin
    scene, layout = layoutscene(resolution=(800, 800))
    ax = layout[1, 1] = LAxis(scene, autolimitaspect=1)

    pie!(ax, 0.1:0.1:1.0, normalize=false)
    scene
end


text:
const KEYWORDS = (
    canvas = BrailleCanvas,
    symbols = ['■'],
    title = "",
    name = "",
    xlabel = "",
    ylabel = "",
    zlabel = "",
    xscale = :identity,
    yscale = :identity,
    labels = true,
    border = :solid,
    width = DEFAULT_WIDTH[],
    height = DEFAULT_HEIGHT[],
    xlim = (0, 0),
    ylim = (0, 0),
    zlim = (0, 0),
    margin = 3,
    padding = 1,
    color = :auto,
    colorbar_lim = (0, 1),
    colorbar_border = :solid,
    colormap = :viridis,
    colorbar = false,
    unicode_exponent = true,
    compact = false,
    blend = true,
    grid = true,
    # internals
    visible = true,
    fix_ar = false,
)

const DESCRIPTION = (
    # NOTE: this named tuple has to stay ordered
    x = "horizontal position for each point",
    y = "vertical position for each point",
    symbols = "characters used to render the bars",
    title = "text displayed on top of the plot",
    name = "current drawing annotation displayed on the right",
    xlabel = "text displayed on the `x` axis of the plot",
    ylabel = "text displayed on the `y` axis of the plot",
    zlabel = "text displayed on the `z` axis (colorbar) of the plot",
    xscale = "`x`-axis scale (`:identity`, `:ln`, `:log2`, `:log10`), or scale function e.g. `x -> log10(x)`",
    yscale = "`y`-axis scale",
    labels = "show plot labels",
    border = "plot bounding box style (`:corners`, `:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, `:none`)",
    margin = "number of empty characters to the left of the whole plot",
    padding = "left and right space between the labels and the canvas",
    color = "choose from (`:green`, `:blue`, `:red`, `:yellow`, `:cyan`, `:magenta`, `:white`, `:normal`, `:auto`), use an integer in `[0-255]`, or provide `3` integers as `RGB` components",
    width = "number of characters per canvas row",
    height = "number of canvas rows",
    xlim = "plotting range for the `x` axis (`(0, 0)` stands for automatic)",
    ylim = "plotting range for the `y` axis",
    zlim = "colormap scaled data range",
    colorbar = "toggle the colorbar",
    colormap = "choose from (`:viridis`, `:cividis`, `:plasma`, `:jet`, `:gray`, `keys(UnicodePlots.COLOR_MAP_DATA)`...), or supply a function `f: (z, zmin, zmax) -> Int(0-255)`, or a vector of RGB tuples",
    colorbar_lim = "colorbar limit",
    colorbar_border = "color bar bounding box style (`:solid`, `:bold`, `:dashed`, `:dotted`, `:ascii`, `:none`)",
    canvas = "type of canvas used for drawing",
    grid = "draws grid-lines at the origin",
    compact = "compact plot labels",
    unicode_exponent = "use `Unicode` symbols for exponents: e.g. `10²⸱¹` instead of `10^2.1`",
    blend = "blend colors on the underlying canvas",
    fix_ar = "fix terminal aspect ratio (experimental)",
    visible = "visible canvas",
)

const Z_DESCRIPTION =
    (:zlabel, :zlim, :colorbar, :colormap, :colorbar_lim, :colorbar_border)

const DEFAULT_KW = (
    # does not have to stay ordered
    :name,
    :title,
    :xlabel,
    :ylabel,
    :zlabel,
    :xscale,
    :yscale,
    :labels,
    :border,
    :width,
    :height,
    :xlim,
    :ylim,
    :zlim,
    :unicode_exponent,
    :compact,
    :blend,
    :grid,
    :padding,
    :margin,
)

const DEFAULT_EXCLUDED = (
    :visible,  # internals
    :fix_ar,  # experimental
    Z_DESCRIPTION...,  # by default for 2D data
)

base_type(x) = replace(string(typeof(x).name.name), "64" => "")

default_with_type(s::Symbol) = (
    if s in keys(KEYWORDS)
        "$s::$(base_type(KEYWORDS[s])) = $(KEYWORDS[s] |> repr)"
    else
        s |> string
    end
)

"""
    keywords([extra]; default = DEFAULT_KW, add = (), exclude = DEFAULT_EXCLUDED, remove = ())

Adds default keywords to a function signature, in a docstring.

# Arguments

- `extra::NamedTuple`: add extra keywords in the form `keyword=value`.
- `default::Tuple`: default `UnicodePlots` keywords.
- `add::Tuple`: add extra symbols, not listed in `default` but present in `KEYWORDS`.
- `remove::Tuple`: remove symbols from `default`.
"""
function keywords(
    extra::NamedTuple = NamedTuple();
    default::Tuple = DEFAULT_KW,
    add::Tuple = (),
    exclude::Tuple = DEFAULT_EXCLUDED,
    remove::Tuple = (),
)
    all_kw = (; KEYWORDS..., extra...)
    candidates = keys(extra) ∪ filter(x -> x ∈ add ∪ default, DEFAULT_KW)
    kw = filter(x -> x ∉ setdiff(exclude ∪ remove, add), candidates)
    @assert allunique(kw)  # extra check
    join((k isa Symbol ? "$k = $(all_kw[k] |> repr)" : k for k in kw), ", ")
end

"""
    arguments([desc]; default = DEFAULT_KW, add = (), exclude = DEFAULT_EXCLUDED, remove = ())

Defines arguments for docstring genreration.

# Arguments

- `desc::NamedTuple`: add argument description in the form `arg=desc`.
- `default::Tuple`: default `UnicodePlots` keywords.
- `add::Tuple`: add extra symbols, not listed in `default` but present in `DESCRIPTION`.
- `remove::Tuple`: remove symbols from `default`.
"""
function arguments(
    desc::NamedTuple = NamedTuple();
    default::Tuple = DEFAULT_KW,
    add::Tuple = (),
    exclude::Tuple = DEFAULT_EXCLUDED,
    remove::Tuple = (),
)
    get_description(n, d) = begin
        str = d[n]
        @assert str[1] === '`' || islowercase(str[1])  # description starts with a lowercase string, or `something`
        @assert str[end] !== '.'  # trailing dot will be added later
        str
    end
    all_desc = (; DESCRIPTION..., desc...)  # desc preempts defaults !
    candidates = keys(desc) ∪ filter(x -> x ∈ add ∪ default, keys(DESCRIPTION))  # desc goes first !
    kw = filter(x -> x ∉ setdiff(exclude ∪ remove, add), candidates)
    @assert allunique(kw)  # extra check
    join(
        ("- **`$(default_with_type(k))`** : $(get_description(k, all_desc))." for k in kw),
        '\n',
    )
end


text:
using Test, Random, DataFrames
using Nanocsv

function rwtest(df; delim=',', na="")
    fname = "nanocsv_test.txt"
    try
        write_csv(fname, df, delim=delim, na=na)
        df2 = read_csv(fname, delim=delim, na=na)
        @test isequal(df, df2)
    finally
        isfile(fname) && rm(fname)
    end
end

function rwthrows(df; err=ArgumentError, delim=',', na="")
    fname = "nanocsv_test.txt"
    try
        write_csv(fname, df, delim=delim, na=na)
        @test_throws err read_csv(fname, delim=delim, na=na)
    finally
        isfile(fname) && rm(fname)
    end
end

@testset "basic read and write" begin
    Random.seed!(1)
    rwtest(DataFrame(rand(Int, 5, 5)))
    rwtest(DataFrame(rand(Int, 5, 5)), delim=';', na="NA")
end

@testset "some more complex test" begin
    df = DataFrame(i = [1,2,missing], f=[missing,1.0, 2.0], m=missing,
                   s1=["s1",",",";"], s2=["s2",missing,missing])
    rwtest(df)
    rwtest(df, delim=';', na="NA")
end

@testset "minimal DataFrames" begin
    rwthrows(DataFrame())
    rwthrows(DataFrame(), delim=';', na="NA")
    rwtest(DataFrame(:a=>[], :b=>[], Symbol("")=>[]))
    rwtest(DataFrame(:a=>[], :b=>[], Symbol("")=>[]), delim=';', na="NA")
    rwtest(DataFrame(Symbol("")=>missing))
    rwtest(DataFrame(Symbol("")=>missing), delim=';', na="NA")
    rwtest(DataFrame(Symbol("")=>["",missing,""]))
    rwtest(DataFrame(Symbol("")=>["",missing,""]), delim=';', na="NA")
end


text:
mutable struct Model{A<:Architecture}
    configuration::ModelConfiguration
    boundary_conditions::ModelBoundaryConditions
    constants::PlanetaryConstants
    eos::EquationOfState
    grid::Grid
    velocities::VelocityFields
    tracers::TracerFields
    pressures::PressureFields
    G::SourceTerms
    Gp::SourceTerms
    forcing  # ::Forcing  # No type so we can set it to nothing while checkpointing.
    stepper_tmp::StepperTemporaryFields
    poisson_solver  # ::PoissonSolver or ::PoissonSolverGPU
    clock::Clock
    output_writers::Array{OutputWriter,1}
    diagnostics::Array{Diagnostic,1}
end

"""
    Model(; kwargs...)

Construct an `Oceananigans.jl` model. The keyword arguments are:

          N (tuple) : model resolution in (x, y, z)
          L (tuple) : model domain in (x, y, z)
         dt (float) : time step
 start_time (float) : start time for the simulation
  iteration (float) : ?
         arch (sym) : architecture (:cpu or :gpu)
  float_type (type) : floating point type for model data (typically Float32 or Float64)
          constants : planetary constants (?)
                eos : equation of state to infer density from temperature and salinity
            forcing : forcing functions for (u, v, w, T, S)
boundary_conditions : boundary conditions
     output_writers : output writer
        diagonstics : diagnostics

        ... and more.
"""
function Model(;
    # Model resolution and domain size
             N,
             L,
    # time step
             Nt,
             Δt,
    # Molecular parameters
             ν = 1.05e-6, νh = ν, νv = ν,
             κ = 1.43e-7, κh = κ, κv = κ,
    # Time stepping
    start_time = 0,
     iteration = 0,
    # Model architecture and floating point precision
          arch = CPU(),
    float_type = Float64,
     constants = Earth(),
    # Equation of State
           eos = LinearEquationOfState(),
    # Forcing and boundary conditions for (u, v, w, T, S)
       forcing = Forcing(nothing, nothing, nothing, nothing, nothing),
    boundary_conditions = ModelBoundaryConditions(),
    # Output and diagonstics
    output_writers = OutputWriter[],
       diagnostics = Diagnostic[]
)

    arch == GPU() && !HAVE_CUDA && throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))

    # Initialize model basics.
    configuration = ModelConfiguration(κh, κv, νh, νv)
             grid = RegularCartesianGrid(float_type, N, L)
            clock = Clock(float_type,start_time,iteration, Nt, Δt)

    # Initialize fields, including source terms and temporary variables.
      velocities = VelocityFields(arch, grid)
         tracers = TracerFields(arch, grid)
       pressures = PressureFields(arch, grid)
               G = SourceTerms(arch, grid)
              Gp = SourceTerms(arch, grid)
     stepper_tmp = StepperTemporaryFields(arch, grid)

     poisson_solver = init_poisson_solver(arch, grid, stepper_tmp.fCC1)

    # Default initial condition
    velocities.u.data .= 0
    velocities.v.data .= 0
    velocities.w.data .= 0
    tracers.S.data .= eos.S₀
    tracers.T.data .= eos.T₀

    Model{typeof(arch)}(configuration, boundary_conditions, constants, eos, grid,
                        velocities, tracers, pressures, G, Gp, forcing,
                        stepper_tmp, poisson_solver, clock, output_writers, diagnostics)
end

arch(model::Model{A}) where A <: Architecture = A
float_type(m::Model) = eltype(model.grid)

add_bcs!(model::Model; kwargs...) = add_bcs(model.boundary_conditions; kwargs...)


text:
#######################################################################
# functions to makle lags and leads
#######################################################################
function F!(df, v, horizons=1, fun=identity)
    # Forward operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("F$(i)$(v)")
        df[:, Fiv] = lead( map(fun, df[!, v]), i)
        push!(newvars, Fiv)
    end
    return newvars
end

function L!(df, v, horizons=1)
    # Lag operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("L$(i)$(v)")
        df[:, Fiv] = lag(df[!, v], i)
        push!(newvars, Fiv)
    end
    return newvars
end

function FΔ!(df, v, horizons=1, fun=identity, baseperiod=0)
    # Forward difference operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("FΔ$(i)$(v)")
        df[:, Fiv] = lead( map(fun, df[!, v]), i) .- lag(map(fun, df[!, v]), baseperiod)
        push!(newvars, Fiv)
    end
    return newvars
end

function LΔ!(df, v, horizons=1, fun=identity)
    # Lag difference operator
    newvars = Symbol[]
    for i in horizons
        Fiv = Symbol("LΔ$(i)$(v)")
        df[:, Fiv] = map(fun, df[!, v]) .-  lag( map(fun, df[!, v]), i)
        push!(newvars, Fiv)
    end
    return newvars
end

function Δ(x, horizons=1, fun=identity)
    return map(fun, x) .-  lag(map(fun, x), horizons)
end


text:
using BulkAndTails, Test

const testfiles = ("pdf_cdf.jl", "fitting.jl")

@testset "BulkAndTails.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
include("fea.jl")

using .FEA: No, Material, Elemento, Estrutura
using Test
using SparseArrays
using LinearAlgebra

@testset verbose = true "Analise" begin
    # =======================================================================
    # Nós
    nos = [No(1, 1, -1; rx=true, ry=true),
           No(2, 3, -1),
           No(3, 5, -1),
           No(4, 7, -1; rx=true, ry=true),
           No(5, 3, 1; fx=10, fy=-5),
           No(6, 5, 1)]

    # Material
    material = Material(1, 100e6)

    # Elementos
    area = 0.0014
    elementos = [Elemento(1, nos[1], nos[2], area, material),
                 Elemento(2, nos[1], nos[5], area, material),
                 Elemento(3, nos[2], nos[3], area, material),
                 Elemento(4, nos[2], nos[6], area, material),
                 Elemento(5, nos[4], nos[3], area, material),
                 Elemento(6, nos[5], nos[2], area, material),
                 Elemento(7, nos[6], nos[3], area, material),
                 Elemento(8, nos[6], nos[4], area, material),
                 Elemento(9, nos[6], nos[5], area, material)]

    # Estrutura
    estrutura = Estrutura(nos, elementos)

    # =======================================================================
    @testset "No" begin
        @test nos[2].id == 2
        @test nos[2].x == 3
        @test nos[2].y == -1
        @test nos[2].fx == 0
        @test nos[2].fy == 0
        @test nos[2].rx == false
        @test nos[2].ry == false

        @test_throws ArgumentError No(0, 3, -1)

        # Métodos
        @test FEA.distancia(nos[1], nos[5]) ≈ 2.8284271247461903

        @test FEA.gls(nos[3]) == [5, 6]

        @test FEA.vetor_forcas(nos[1]) == [0, 0]
        @test FEA.vetor_forcas(nos[5]) == [10, -5]

        @test FEA.vetor_apoios(nos[1]) == [true, true]
        @test FEA.vetor_apoios(nos[2]) == [false, false]

        nos_col = [No(1, 1, 1), No(2, 2, 2), No(3, 3, 3)]
        @test FEA.sao_colineares(nos_col...) == true
        nos_col[1].x = 1.0001
        @test FEA.sao_colineares(nos_col...) == false

        @test nos_col[1].fx == 0
        @test nos_col[1].fy == 0
        @test nos_col[1].rx == false
        @test nos_col[1].ry == false
        FEA.carregar_apoiar_no!(nos_col, 1.1, 1.1; fx=5, fy=6, rx=true, ry=true)
        @test nos_col[1].fx == 5
        @test nos_col[1].fy == 6
        @test nos_col[1].rx == true
        @test nos_col[1].ry == true
    end

    # =======================================================================
    @testset "Material" begin
        @test material.id == 1
        @test material.E == 100e6

        @test_throws ArgumentError Material(-1, 100e6)
        @test_throws ArgumentError Material(1, -100e6)
    end

    # =======================================================================
    @testset "Elementos" begin
        @test elementos[2].id == 2
        @test elementos[2].no₁ == nos[1]
        @test elementos[2].no₂ == nos[5]
        @test elementos[2].area ≈ area
        @test elementos[2].material == material

        # Métodos
        @test FEA.comprimento(elementos[2]) ≈ 2.8284271247461903

        @test FEA.gls(elementos[2]) == [1, 2, 9, 10]

        @test FEA.vetor_forcas(elementos[2]) == [0, 0, 10, -5]
        @test FEA.vetor_forcas(elementos[9]) == [0, 0, 10, -5]
        @test FEA.vetor_forcas(elementos[3]) == [0, 0, 0, 0]

        @test FEA.angulo(elementos[4]) ≈ pi / 4
        @test FEA.angulo(elementos[9]) ≈ pi
        @test FEA.angulo(elementos[8]) ≈ -pi / 4

        @test FEA.matriz_rotacao(elementos[4]) ≈ [√2/2 √2/2  0       0
                                                      -√2/2 √2/2  0       0
                                                       0       0  √2/2 √2/2
                                                       0       0 -√2/2 √2/2]

        @test FEA.k_local(elementos[4]) ≈ [49497.474683058324  0 -49497.474683058324   0
                                                           0  0                   0   0
                                         -49497.474683058324  0  49497.474683058324   0
                                                           0  0                   0   0]

        @test FEA.k(elementos[4]) ≈ [24748.737341529166 24748.737341529162 -24748.737341529166 -24748.737341529162
                                                 24748.73734152916 24748.737341529155 -24748.73734152916 -24748.737341529155
                                                -24748.737341529166 -24748.737341529162 24748.737341529166 24748.737341529162
                                                -24748.73734152916 -24748.737341529155 24748.73734152916 24748.737341529155]

        # Elementos sobreposição
        nos_sob = [No(1,0,0),
                   No(2,50,50),
                   No(3,10,10),
                   No(4,30,30),
                   No(5,-10,-10),
                   No(6,60,60),
                   No(7,-20,-20),
                   No(8,40,30),
                   No(9,10,40)]

        els_sob = [Elemento(1, nos_sob[1], nos_sob[2], 1, material),
                   Elemento(2, nos_sob[3], nos_sob[4], 1, material),
                   Elemento(3, nos_sob[3], nos_sob[5], 1, material),
                   Elemento(4, nos_sob[5], nos_sob[6], 1, material),
                   Elemento(5, nos_sob[2], nos_sob[1], 1, material),
                   Elemento(6, nos_sob[5], nos_sob[7], 1, material),
                   Elemento(7, nos_sob[8], nos_sob[9], 1, material)]

        @test FEA.tem_sobreposicao(els_sob[1], els_sob[2]) == true
        @test FEA.tem_sobreposicao(els_sob[1], els_sob[3]) == true
        @test FEA.tem_sobreposicao(els_sob[1], els_sob[4]) == true
        @test FEA.tem_sobreposicao(els_sob[1], els_sob[5]) == true
        @test FEA.tem_sobreposicao(els_sob[1], els_sob[6]) == false
        @test FEA.tem_sobreposicao(els_sob[1], els_sob[7]) == false

    end

    # =======================================================================
    @testset "Estrutura" begin
        @test estrutura.nos == nos
        @test estrutura.elementos == elementos

        # Métodos
        @test FEA.num_gls_estrutura(estrutura) == 12

        @test FEA.gls_livres(estrutura) == [3, 4, 5, 6, 9, 10, 11, 12]

        @test FEA.vetor_forcas(estrutura, true) == [0, 0, 0, 0, 0, 0, 0, 0, 10, -5, 0, 0]
        @test FEA.vetor_forcas(estrutura, false) == [0, 0, 0, 0, 10, -5, 0, 0]

        @test FEA.vetor_apoios(estrutura) == [1, 2, 7, 8]

        kt = sparse([164748.73734152917 24748.73734152916 -70000.0 0.0 -2.624579619658251e-28 4.2862637970157364e-12 -24748.737341529166 -24748.737341529162
                    24748.737341529155 94748.73734152915 0.0 0.0 4.2862637970157364e-12 -70000.0 -24748.73734152916 -24748.737341529155
                    -70000.0 0.0 140000.0 -4.2862637970157364e-12 0.0 0.0 -2.624579619658251e-28 4.2862637970157364e-12
                    0.0 0.0 -4.2862637970157364e-12 70000.0 0.0 0.0 4.2862637970157364e-12 -70000.0
                    -2.624579619658251e-28 4.2862637970157364e-12 0.0 0.0 94748.73734152917 24748.73734152916 -70000.0 0.0
                    4.2862637970157364e-12 -70000.0 0.0 0.0 24748.737341529155 94748.73734152915 0.0 0.0
                    -24748.737341529166 -24748.737341529162 -2.624579619658251e-28 4.2862637970157364e-12 -70000.0 0.0 119497.47468305833 -3.637978807091713e-12
                    -24748.73734152916 -24748.737341529155 4.2862637970157364e-12 -70000.0 0.0 0.0 -3.637978807091713e-12 119497.4746830583])

        @test FEA.k(estrutura) ≈ kt

        @test FEA.deslocamentos(estrutura) ≈ [0.0, 0.0, 4.761904761904763e-5, -0.00019817906943235845, 2.3809523809523814e-5, -7.528001090665546e-5, 0.0, 0.0, 0.0002696076408609297, -0.0002696076408609298, 0.00012675049800378677, -7.528001090665545e-5]

        @test FEA.forcas_internas(estrutura) ≈ [3.3340206284604434, 0.0012720448909134228, -1.6664023403833994, 7.0710761922049326, -1.6665674915761208, -4.999137587831104, 0.00018544477509152557, -7.07105622101196, -9.999339096309345]
    end

end

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


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

using Markdown
using InteractiveUtils

# ╔═╡ 02c81da6-3975-11ec-08d6-a35a7c2772d6
md"""
# The `lib/` folder

Genie makes it very easy to automatically load your existing Julia code (modules, files, etc) into the app, outside of the standard Genie MVC structure. You simply need to add your files and folders into the `lib/` folder.

---
**HEADS UP**

* If the `lib/` folder does not exist, just create it yourself: `julia> mkdir("lib")`
* Genie does not `include` the files placed within the `lib/` folder, but _recursively_ adds all the folders to the `LOAD_PATH` - so you can `include` the files yourself, as needed -- Genie/Julia will know were to find them in the LOAD_PATH.
* Files within `lib/` are not added to the `Revise` queue so they are not automatically reloaded by Genie if changed. If you make changes/add/remove files in `lib/` you need to restart the app or manually add them to be watched and reloaded by `Revise`. The reason being that the `lib/` folder is meant for legacy code. If you will actively develop the code, the recommended way is to use Genie's MVC structure, adding Models, Controllers, and Views as needed.

---"""

# ╔═╡ Cell order:
# ╟─02c81da6-3975-11ec-08d6-a35a7c2772d6


text:
############################Chronologies For FeedForward###################################
@doc raw"""
    Synchronize(periods::Int)
Defines the co-ordination of time between Two stages.

# Arguments
- `periods::Int`: Number of time periods to grab data from
"""
struct Synchronize <: FeedForwardChronology
    periods::Int
    function Synchronize(; periods)
        new(periods)
    end
end
# TODO: Add DocString
"""
    RecedingHorizon(period::Int)
"""
struct RecedingHorizon <: FeedForwardChronology
    periods::Int
    function RecedingHorizon(; periods::Int = 1)
        new(periods)
    end
end

struct Consecutive <: FeedForwardChronology end

function check_chronology!(sim::Simulation, key::Pair, sync::Synchronize)
    from_stage = get_stage(sim, key.first)
    to_stage = get_stage(sim, key.second)
    from_stage_horizon = sim.sequence.horizons[key.first]
    to_stage_horizon = sim.sequence.horizons[key.second]
    from_stage_interval = get_stage_interval(sim, key.first)
    to_stage_interval = get_stage_interval(sim, key.second)

    from_stage_resolution =
        IS.time_period_conversion(PSY.get_forecasts_resolution(from_stage.sys))
    @debug from_stage_resolution, to_stage_interval
    # How many times the second stages executes per solution retireved from the from_stage.
    # E.g. from_stage_resolution = 1 Hour, to_stage_interval = 5 minutes => 12 executions per solution
    to_stage_executions_per_solution = Int(from_stage_resolution / to_stage_interval)
    # Number of periods in the horizon that will be synchronized between the from_stage and the to_stage
    from_stage_sync = sync.periods

    if from_stage_sync > from_stage_horizon
        throw(IS.ConflictingInputsError("The lookahead length $(from_stage_horizon) in stage is insufficient to syncronize with $(from_stage_sync) feedforward periods"))
    end

    if (from_stage_sync % to_stage_executions_per_solution) != 0
        throw(IS.ConflictingInputsError("The current configuration implies $(from_stage_sync / to_stage_executions_per_solution) executions of $(key.second) per execution of $(key.first). The number of Synchronize periods $(sync.periods) in stage $(key.first) needs to be a mutiple of the number of stage $(key.second) execution for every stage $(key.first) interval."))
    end

    return
end

function check_chronology!(sim::Simulation, key::Pair, ::Consecutive)
    from_stage_horizon = sim.sequence.horizons[key.first]
    from_stage_interval = get_stage_interval(sim, key.first)
    if from_stage_horizon != from_stage_interval
        @warn("Consecutive Chronology Requires the same interval and horizon, the parameter horizon = $(from_stage_horizon) in stage $(key.first) will be replaced with $(from_stage_interval). If this is not the desired behviour consider changing your chronology to RecedingHorizon")
    end
    sim.sequence.horizons[key.first] = get_stage_interval(sim, key.first)
    return
end

check_chronology!(sim::Simulation, key::Pair, ::RecedingHorizon) = nothing

function check_chronology!(
    sim::Simulation,
    key::Pair,
    ::T,
) where {T <: FeedForwardChronology}
    error("Chronology $(T) not implemented")
    return
end

############################FeedForward Definitions########################################

struct UpperBoundFF <: AbstractAffectFeedForward
    variable_from_stage::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function UpperBoundFF(
        variable_from_stage::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(variable_from_stage), Symbol.(affected_variables), cache)
    end
end

function UpperBoundFF(; variable_from_stage, affected_variables)
    return UpperBoundFF(variable_from_stage, affected_variables, nothing)
end

get_variable_from_stage(p::UpperBoundFF) = p.binary_from_stage

struct RangeFF <: AbstractAffectFeedForward
    variable_from_stage_ub::Symbol
    variable_from_stage_lb::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function RangeFF(
        variable_from_stage_ub::AbstractString,
        variable_from_stage_lb::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(
            Symbol(variable_from_stage_ub),
            Symbol(variable_from_stage_ub),
            Symbol.(affected_variables),
            cache,
        )
    end
end

function RangeFF(; variable_from_stage_ub, variable_from_stage_lb, affected_variables)
    return RangeFF(
        variable_from_stage_ub,
        variable_from_stage_lb,
        affected_variables,
        nothing,
    )
end

get_bounds_from_stage(p::RangeFF) = (p.variable_from_stage_lb, p.variable_from_stage_lb)

struct SemiContinuousFF <: AbstractAffectFeedForward
    binary_from_stage::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function SemiContinuousFF(
        binary_from_stage::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(binary_from_stage), Symbol.(affected_variables), cache)
    end
end

function SemiContinuousFF(; binary_from_stage, affected_variables)
    return SemiContinuousFF(binary_from_stage, affected_variables, nothing)
end

get_binary_from_stage(p::SemiContinuousFF) = p.binary_from_stage
get_affected_variables(p::AbstractAffectFeedForward) = p.affected_variables

struct IntegralLimitFF <: AbstractAffectFeedForward
    variable_from_stage::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function IntegralLimitFF(
        variable_from_stage::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(variable_from_stage), Symbol.(affected_variables), cache)
    end
end

function IntegralLimitFF(; variable_from_stage, affected_variables)
    return IntegralLimitFF(variable_from_stage, affected_variables, nothing)
end

get_variable_from_stage(p::IntegralLimitFF) = p.variable_from_stage

####################### Feed Forward Affects ###############################################

@doc raw"""
        ub_ff(psi_container::PSIContainer,
              cons_name::Symbol,
              param_reference::UpdateRef,
              var_name::Symbol)

Constructs a parametrized upper bound constraint to implement feedforward from other models.
The Parameters are initialized using the uppper boundary values of the provided variables.

# Constraints
``` variable[var_name, t] <= param_reference[var_name] ```

# LaTeX

`` x \leq param^{max}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the PJ.ParameterRef used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function ub_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    variable = get_variable(psi_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]

    @assert axes[2] == time_steps
    container = add_param_container!(psi_container, param_reference, set_name)
    param_ub = get_parameter_array(container)
    con_ub = add_cons_container!(psi_container, ub_name, set_name, time_steps)

    for name in axes[1]
        value = JuMP.upper_bound(variable[name, 1])
        param_ub[name] = PJ.add_parameter(psi_container.JuMPmodel, value)
        for t in axes[2]
            con_ub[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                variable[name, t] <= param_ub[name]
            )
        end
    end

    return
end

@doc raw"""
        range_ff(psi_container::PSIContainer,
                        cons_name::Symbol,
                        param_reference::NTuple{2, UpdateRef},
                        var_name::Symbol)

Constructs min/max range parametrized constraint from device variable to include feedforward.

# Constraints

``` param_reference[1][var_name] <= variable[var_name, t] ```
``` variable[var_name, t] <= param_reference[2][var_name] ```

where r in range_data.

# LaTeX

`` param^{min} \leq x ``
`` x \leq param^{max}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* param_reference::NTuple{2, UpdateRef} : Tuple with the lower bound and upper bound parameter reference
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
"""
function range_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    param_reference::NTuple{2, UpdateRef},
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")

    variable = get_variable(psi_container, var_name)
    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps

    #Create containers for the constraints
    container_lb = add_param_container!(psi_container, param_reference[1], set_name)
    param_lb = get_parameter_array(container_lb)
    container_ub = add_param_container!(psi_container, param_reference[2], set_name)
    param_ub = get_parameter_array(container_ub)

    #Create containers for the parameters
    con_lb = add_cons_container!(psi_container, lb_name, set_name, time_steps)
    con_ub = add_cons_container!(psi_container, ub_name, set_name, time_steps)

    for name in axes[1]
        param_lb[name] =
            PJ.add_parameter(psi_container.JuMPmodel, JuMP.lower_bound(variable[name, 1]))
        param_ub[name] =
            PJ.add_parameter(psi_container.JuMPmodel, JuMP.upper_bound(variable[name, 1]))
        for t in axes[2]
            con_ub[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                variable[name, t] <= param_ub[name]
            )
            con_lb[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                variable[name, t] >= param_lb[name]
            )
        end
    end

    return
end

@doc raw"""
            semicontinuousrange_ff(psi_container::PSIContainer,
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    param_reference::UpdateRef)

Constructs min/max range constraint from device variable with parameter setting.

# Constraints
If device min = 0:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

Otherwise:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

``` variable[var_name, t] >= r[2].min*param_reference[var_name] ```

where r in range_data.

# LaTeX

`` 0.0 \leq x^{var} \leq r^{max} x^{param}, \text{ for } r^{min} = 0 ``

`` r^{min} x^{param} \leq x^{var} \leq r^{min} x^{param}, \text{ otherwise } ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* param_reference::UpdateRef : UpdateRef of the parameter
"""
function semicontinuousrange_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")

    variable = get_variable(psi_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps
    container = add_param_container!(psi_container, param_reference, set_name)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    con_ub = add_cons_container!(psi_container, ub_name, set_name, time_steps)
    con_lb = add_cons_container!(psi_container, lb_name, set_name, time_steps)

    for name in axes[1]
        ub_value = JuMP.upper_bound(variable[name, 1])
        lb_value = JuMP.lower_bound(variable[name, 1])
        param[name] = PJ.add_parameter(psi_container.JuMPmodel, 1.0)
        for t in axes[2]
            con_ub[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                variable[name, t] <= ub_value * param[name]
            )
            con_lb[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                variable[name, t] >= lb_value * param[name]
            )
        end
    end

    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            JuMP.set_lower_bound(v, 0.0)
        end
    end

    return
end

@doc raw"""
        integral_limit_ff(psi_container::PSIContainer,
                        cons_name::Symbol,
                        param_reference::UpdateRef,
                        var_name::Symbol)

Constructs a parametrized integral limit constraint to implement feedforward from other models.
The Parameters are initialized using the upper boundary values of the provided variables.

# Constraints
``` sum(variable[var_name, t] for t in time_steps)/length(time_steps) <= param_reference[var_name] ```

# LaTeX

`` \sum_{t} x \leq param^{max}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the PJ.ParameterRef used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function integral_limit_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "integral_limit")
    variable = get_variable(psi_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]

    @assert axes[2] == time_steps
    container_ub = add_param_container!(psi_container, param_reference, set_name)
    param_ub = get_parameter_array(container_ub)
    con_ub = add_cons_container!(psi_container, ub_name, set_name)

    for name in axes[1]
        value = JuMP.upper_bound(variable[name, 1])

        param_ub[name] = PJ.add_parameter(psi_container.JuMPmodel, value)
        con_ub[name] = JuMP.@constraint(
            psi_container.JuMPmodel,
            sum(variable[name, t] for t in time_steps) / length(time_steps) <=
            param_ub[name]
        )
    end
end

########################## FeedForward Constraints #########################################
function feedforward!(
    psi_container::PSIContainer,
    device_type::Type{T},
    ff_model::Nothing,
) where {T <: PSY.Component}
    return
end

function feedforward!(
    psi_container::PSIContainer,
    device_type::Type{I},
    ff_model::UpperBoundFF,
) where {I <: PSY.StaticInjection}
    for prefix in get_affected_variables(ff_model)
        var_name = variable_name(prefix, I)
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        ub_ff(psi_container, constraint_name(feedforward, I), parameter_ref, var_name)
    end
end

function feedforward!(
    psi_container::PSIContainer,
    ::Type{T},
    ff_model::SemiContinuousFF,
) where {T <: PSY.StaticInjection}
    bin_var = variable_name(get_binary_from_stage(ff_model), T)
    parameter_ref = UpdateRef{JuMP.VariableRef}(bin_var)
    for prefix in get_affected_variables(ff_model)
        var_name = variable_name(prefix, T)
        semicontinuousrange_ff(
            psi_container,
            constraint_name(FEEDFORWARD_BIN, T),
            parameter_ref,
            var_name,
        )
    end
end

function feedforward!(
    psi_container::PSIContainer,
    ::Type{T},
    ff_model::IntegralLimitFF,
) where {T <: PSY.StaticInjection}
    for prefix in get_affected_variables(ff_model)
        var_name = variable_name(prefix, T)
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        integral_limit_ff(
            psi_container,
            constraint_name(FEEDFORWARD, T),
            parameter_ref,
            var_name,
        )
    end
end

#########################FeedForward Variables Updating#####################################
# This makes the choice in which variable to get from the results.
function get_stage_variable(
    ::RecedingHorizon,
    stages::Pair{Stage{T}, Stage{T}},
    device_name::AbstractString,
    var_ref::UpdateRef,
) where {T <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    step = axes(variable)[2][1]
    return JuMP.value(variable[device_name, step])
end

function get_stage_variable(
    ::Consecutive,
    stages::Pair{Stage{T}, Stage{V}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, V <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    step = axes(variable)[2][get_end_of_interval_step(stages.first)]
    return JuMP.value(variable[device_name, step])
end

function get_stage_variable(
    ::Synchronize,
    stages::Pair{Stage{T}, Stage{T}},
    device_name::String,
    var_ref::UpdateRef,
) where {T <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    step = axes(variable)[2][stages.second.internal.execution_count + 1]
    return JuMP.value(variable[device_name, step])
end

function feedforward_update(
    sync::FeedForwardChronology,
    param_reference::UpdateRef{JuMP.VariableRef},
    param_array::JuMPParamArray,
    to_stage::Stage,
    from_stage::Stage,
)
    for device_name in axes(param_array)[1]
        var_value =
            get_stage_variable(sync, (from_stage => to_stage), device_name, param_reference)
        PJ.fix(param_array[device_name], var_value)
    end
end


text:
"""
    contract(A, idxA, B, idxB)

contract the symmetric tensor A, B. The to-be-contracted
indexes are shown with negative integers while the remaining
indexes for the result are shown with positive integers. For
example if there are l negative numbers in both A and B (should
be numbered from -1 to -l) then the rest of the indexes in A and B
combined should be 1:N+M-2n.
"""
function contract(A     :: AbstractSymTensor{S,T1,N},
                  idxA  :: NTuple{N, Int},
                  B     :: AbstractSymTensor{S,T2,M},
                  idxB  :: NTuple{M, Int}) where{S,T1,T2, N, M}

    remsA, consA, tofinalsA = _contract_index_perm(idxA)
    remsB, consB, tofinalsB = _contract_index_perm(idxB)

    ## TODO: techniqually here the contraction vector space should be made
    ## once! But now is made twice, so can we do better!

    # TODO: check
    # compatibility println(remsA, consA, remsB, consB)
    if length(remsA) == 0
        _A = SymVector(A, consA, false)
        if length(remsB) == 0
            _B = SymVector(B, consB, true)
            return _A * _B
        else
            _B = SymMatrix(B, consB, remsB)
            return permutelegs(splitleg(_A*_B, 2, B.legs[remsB]),
                               invperm(tofinalsB))
        end
    end

    _A = SymMatrix(A, remsA, consA)

    if length(remsB) == 0
        _B = SymVector(B, consB, true)
        return permutelegs(splitleg(_A*_B, 1, A.legs[remsA]),
                           invperm(tofinalsA))
    end

    _B = SymMatrix(B, consB, remsB)
    return permutelegs(SymTensor(_A * _B, A.space[remsA], B.space[remsB]),
                       invperm([tofinalsA; tofinalsB]))
end

function contract(A     :: AbstractArray{T1,N},
                  idxA  :: NTuple{N, Int},
                  B     :: AbstractArray{T2,M},
                  idxB  :: NTuple{M, Int}) where{S,T1,T2, N, M}
    remsA, consA, tofinalsA = _contract_index_perm(idxA)
    remsB, consB, tofinalsB = _contract_index_perm(idxB)
    permA = [remsA; consA]
    permB = [consB; remsB]

    if length(remsA) == 0
        _A = reshape(permutedims(A, permA), prod(size(A)[permA]))
        if length(remsB) == 0
            _B = reshape(permutedims(B, permB), prod(size(B)[permB]))
            return sum(_A .* _B)
        else
            _B = reshape(permutedims(B, permB), prod(size(B)[consB]), prod(size(B)[remsB]))
            return permutedims(reshape(transpose(_A)*_B, size(B)[remsB]...),
                               invperm(tofinalsB))
        end
    end

    _A = reshape(permutedims(A, permA), prod(size(A)[remsA]), prod(size(A)[consA]))
    if length(remsB) == 0
        _B = reshape(permutedims(B, permB), prod(size(B)[permB]))
        return permutedims(reshape(_A*_B, size(A)[remsA]...),
                           invperm(tofinalsA))
    end
    _B = reshape(permutedims(B, permB), prod(size(B)[consB]), prod(size(B)[remsB]))
    permutedims(reshape(_A*_B, size(A)[remsA]..., size(B)[remsB]...),
                invperm([tofinalsA; tofinalsB]))
end

# auxillary function to find indexes for contract
function _contract_index_perm(indexset; mode::Symbol=:MINUS)
    rems = Int[]
    cons = Int[]
    cons_order = Int[]
    tofinals = Int[]
    for i in eachindex(indexset)
        idx = indexset[i]
        if idx < 0
            push!(cons, i)
            push!(cons_order, idx)
        else
            push!(rems, i)
            push!(tofinals, idx)
        end
    end
    perm = sortperm(tofinals)
    rems[perm], cons[sortperm(cons_order, by=x->abs(x))], tofinals[perm]
end

function SymVector(A::AbstractSymTensor,
                   perm::Vector{Int},
                   _dual::Bool)
    N = rank(A)
    sort(perm) == collect(1:N) ||
        error("Incorrect index set for conversion to SymVector!")

    pA = permutelegs(A, perm)
    fuselegs(pA, 1, N, _dual)
end

function SymMatrix(A::AbstractSymTensor,
                   rowidxs::Vector{Int},
                   colidxs::Vector{Int})

    N = rank(A)
    idxperm = [rowidxs; colidxs]
    sort(idxperm) == collect(1:N) ||
        error("Incorrect index set for conversion to SymMatrix!")
    length(rowidxs) == 0 && error("empty row for matrix not allowed!")
    length(colidxs) == 0 && error("empty col for matrix not allowed!")
    n = length(rowidxs)

    T = eltype(A)
    S = vtype(A)

    rowdualinfo = isdual.(A.space[rowidxs])
    coldualinfo = isdual.(A.space[colidxs])
    csects = [Sector(x[idxperm]) for x in sectors(A)]
    sperm1 = sortperm(csects)
    fsects = Vector{Sector{S, 2}}(undef, length(csects))
    for i in eachindex(csects)
        sector = csects[sperm1][i]
        c1 = sum(Sector(sector[1:n]), rowdualinfo)
        c2 = inv(sum(Sector(sector[n+1:N]), coldualinfo))
        fsects[i] = Sector(c1, c2)
    end
    sperm2 = sortperm(fsects)
    semits = collect(onlysemitokens(A.blocks))[sperm1][sperm2]

    space = (fuse(false, A.space[rowidxs]), fuse(true, A.space[colidxs]))

    sects, sizes = _allsectorsandsizes(A.charge, space)

    blocks = SortedDict{Sector{S, 2}, Matrix{T}}()
    pointer = 1
    for index in 1:length(sects)
        blk = Matrix{T}(undef, sizes[index])

        rlim, clim = sizes[index]
        pc = 1
        while pc <= clim
            pr = 1
            fdc = 0
            while pr <= rlim
                Ablk = permutedims(A.blocks[semits[pointer]], idxperm)
                s = size(Ablk)
                fdr = prod(s[1:n])
                fdc = prod(s[n+1:N])
                blk[pr:pr+fdr-1, pc:pc+fdc-1] =
                    reshape(Ablk, fdr, fdc)
                pr += fdr
                pointer += 1
            end
            pc += fdc
        end
        sector = sects[index]
        blocks[sector] = blk
    end
    SymMatrix{S,T}(A.charge, space, blocks)
end

function SymMatrix(A::Array,
                   rowidxs::Vector{Int},
                   colidxs::Vector{Int})

    N = ndims(A)
    idxperm = [rowidxs; colidxs]
    sort(idxperm) == collect(1:N) ||
        error("Incorrect index set for conversion to SymMatrix!")
    length(rowidxs) == 0 && error("empty row for matrix not allowed!")
    length(colidxs) == 0 && error("empty col for matrix not allowed!")
    T = eltype(A)
    reshape(permutedims(A, idxperm), prod(size(A)[rowidxs]), prod(size(A)[colidxs]))
end

function SymTensor(A     :: SymMatrix,
                   rlegs :: NTuple{N, VectorSpace{S}},
                   clegs :: NTuple{M, VectorSpace{S}}) where {S, N, M}
    T = eltype(A)

    blocks = SortedDict{Sector{S, N+M}, Array{T, N+M}}()
    csects = sectors(A)

    rpats, rsizes = Vector{Sector{S, M}}(), Vector{NTuple{M, Int}}()
    cpats, csizes = Vector{Sector{S, M}}(), Vector{NTuple{M, Int}}()
    isduals = isdual.(A.space)
    for (s,blk) in A.blocks
        c1, c2 = s.charges
        c1 = isdual(A.space[1]) ? inv(c1) : c1
        c2 = isdual(A.space[2]) ? inv(c2) : c2
        rpats, rsizes = _allsectorsandsizes(c1, rlegs)
        cpats, csizes = _allsectorsandsizes(c2, clegs)

        pc = 0
        for cpatidx in eachindex(cpats)
            csl = prod(csizes[cpatidx])
            pr = 0
            for rpatidx in eachindex(rpats)
                rsl = prod(rsizes[rpatidx])
                blocks[Sector(rpats[rpatidx]..., cpats[cpatidx]...)] =
                    reshape(blk[pr+1:pr+rsl, pc+1:pc+csl],
                            rsizes[rpatidx]..., csizes[cpatidx]...)
                pr += rsl
            end
            pc += csl
        end
    end

    #TODO: now check to see if the new legs can fuse into the original leg
    space = (rlegs...,clegs...)

    SymTensor{S,T,N+M}(A.charge, space, blocks)
end

function *(A::SymVector, B::SymVector)
    isdual(A.space[1], B.space[1]) ||
        error("not contractible!", A.space[1], " and ", B.space[1])

    semits = collect(onlysemitokens(A.blocks))
    return sum(A.blocks[semits[1]] .* B.blocks[semits[1]])
end

function *(A::AbstractSymMatrix,
           B::AbstractSymMatrix)

    isdual(A.space[2], B.space[1]) ||
        error("not contractible! ", A.space[2], " and ", B.space[1])

    S = vtype(A)
    T = promote_type(eltype(A), eltype(B))
    T == Union{} &&
        error("can't promote_type! $(eltype(A)) , $(eltype(B))")
    ## NOTE: assume the set of charges for the vector space to be
    ## contracted is {c1,..,ci,...,cn} then each sector in A is given
    ## by (CA+ci, -ci) and each sector in B by (ci, -ci+CB) and since
    ## both are sorted based on the last charge in sector, then they
    ## are sorted the same. So we simply need to multiply them here
    ## but take into account that it is possible for the matrices to
    ## have different number of sectors allowed due to the other
    ## vector space they map from/to.

    ## NOTE: After contraction it is possible that some new sectors
    ## (that weren't allowed with the contraction indexes) are now
    ## allowed. So they have to be manually added and initialized to
    ## zero. I am actually not very happy with this; originally I
    ## decided to not include these new sectors at all but then
    ## realized at some instances the existance of them is necessary,
    ## well they are allowed sectors anyways, right! Still don't know
    ## what is the best approach.

    ## TODO: Also it would be nice if there is the option of dropping
    ## the charges that only show up in complete zero sectors! This is
    ## particularly useful is Gutzwiller projection thing

    charge = A.charge + B.charge
    space = (A.space[1], B.space[2])

    sects, sizes = _allsectorsandsizes(charge, space)
    n_sectors = length(sects)
    blocks = SortedDict{Sector{S, 2}, Matrix{T}}()

    semitsA = collect(onlysemitokens(A.blocks))
    semitsB = collect(onlysemitokens(B.blocks))
    bmax = length(semitsB)
    amax = length(semitsA)
    a, b = 1, 1

    for i=1:n_sectors
        c1, c2 = sects[i]
        _c1 = isdual(A.space[1]) ? inv(c1) : c1
        _c2 = isdual(B.space[2]) ? inv(c2) : c2
        if isdual(A.space[2])
            sectorA = Sector(c1, inv(A.charge)+_c1)
            sectorB = Sector(B.charge-_c2, c2)
        else
            sectorA = Sector(c1, A.charge-_c1)
            sectorB = Sector(inv(B.charge)+_c2, c2)
        end
        if haskey(A.blocks, sectorA) && haskey(B.blocks, sectorB)
            blocks[sects[i]] = A.blocks[sectorA] * B.blocks[sectorB]
        else
            blocks[sects[i]] = zeros(T, sizes[i])
        end
    end
    #SymMatrix(charge, space, blocks)
    SymMatrix{S,T}(charge, space, blocks)
end

###TODO: This is the above function with a version of fuselegs (that fuses both of them at once!)
function SymMatrix2(A::AbstractSymTensor{T, N},
                    rowidxs::Vector{Int},
                    colidxs::Vector{Int}) where {T<:Number, N}

    idxperm = [rowidxs; colidxs]
    sort(idxperm) == collect(1:N) || error("Can't convert to SymMatrix!")
    length(rowidxs) == 0 && error("Zero row for matrix not allowed!")
    length(colidxs) == 0 && error("Zero col for matrix not allowed!")

    pA = permutelegs(A, idxperm)

    fusedsectors = NTuple{Int, Int}[]
    patranges = NTuple{UnitRange{Int}, UnitRange{Int}}[]

    #     csects = sten.sects
    #     fchrdict = Dict{Int, FusedCharge{n}}()
    #     signs = Tuple(sten.legs[c].sign for c in l:l+n-1)
    #     for i in eachindex(csects)
    #         sect = [csects[i]...]

    #         pat = Tuple(sect[l:l+n-1])
    #         spat = sign .* signs .* pat
    #         fchr = sign * sum(signs .* [sect[c] for c in l:l+n-1])
    #         patrange = UnitRange{Int}(1,1)
    #         if haskey(fchrdict, fchr)
    #             if haskey(fchrdict[fchr].pats, spat)
    #                 patrange = fchrdict[fchr].pats[spat]
    #             else
    #                 fdim = prod([getdim(sten.legs[c+l-1], pat[c]) for c in 1:n])
    #                 dim = fchrdict[fchr].dim
    #                 fchrdict[fchr].dim += fdim
    #                 patrange = dim+1:dim+fdim
    #                 fchrdict[fchr].pats[spat] = patrange
    #             end
    #         else
    #             fdim = prod([getdim(sten.legs[c+l-1], pat[c]) for c in 1:n])
    #             patrange = 1:fdim
    #             fchrdict[fchr] = FusedCharge(fchr, fdim, Dict(spat => patrange))
    #         end

    #         fsect = vcat( sect[1:l-1], fchr, sect[l+n:end])
    #         push!(fusedsectors, Tuple(fsect))
    #         push!(patranges, patrange)
    #     end

    #     fusedsectperm = sortperm(fusedsectors)
    #     new_sectors, refs = uniquesorted(fusedsectors[fusedsectperm])
    #     refs = refs[invperm(fusedsectperm)]

    #     # the infos are: new_sectors, refs, patranges

    #     # make the new leg
    #     fleg = U1Space(sign, unzip([(k, fchrdict[k].dim) for k in sort(collect(keys(fchrdict)))])...)

    #     new_legs = Tuple(vcat([sten.legs[1:l-1]...], fleg, [sten.legs[l+n:end]...]))

    #     new_nzblks = Array{T, N-n+1}[]
    #     for sect in new_sectors
    #         push!(new_nzblks,Array{T, N-n+1}(undef, [getdim(new_legs[c], sect[c]) for c in 1:N-n+1]...))
    #     end
    #     #indexing = [1:end for i=1:N-n+1]...
    #     for i in eachindex(sten.nzblks)
    #         s = size(sten.nzblks[i])
    #         new_nzblks[refs[i]][[1:s[i] for i=1:l-1]..., patranges[i],[1:s[i] for i=l+n:N]...] =
    #             reshape(sten.nzblks[i], s[1:l-1]..., prod(s[l:l+n-1]), s[l+n:end]...)
    #     end

    #     SymTensor(sten.charge, new_legs, new_sectors, new_nzblks)

    fuselegs(
        fuselegs(psten, -1, length(rowidxs)+1, length(colidxs)),
        +1, 1, length(rowidxs))

end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
using JSON
Nproc=4

rmprocs(workers())
addprocs(Nproc)

println("workers: ",length(workers()))

@everywhere using Base.Test
@everywhere import Base.convert
@everywhere using MAT
@everywhere using ForwardDiff
@everywhere using Polynomials
@everywhere import ForwardDiff.DiffBase
@everywhere using Optim
@everywhere gc(true)

# local
@everywhere path="/home/emma/github/fit_behaviour_julia/"
# cluster
# @everywhere path="/home/hcli64/hcli64751/rats/fit_behaviour_julia/"
#
@everywhere include(path*"/functions/functions_julia.jl")
@everywhere include(path*"/functions/parameter_functions.jl")
@everywhere include(path*"/functions/utils.jl")

rat_num = 5
path_load=path*"scripts_rats/"
filename="data/ordered_processed_dataset2_1sec.json"

data = Dict()
open(path_load*filename, "r") do f
    global data
    dicttxt = readstring(f)  # file information to string
    data=JSON.parse(dicttxt)  # parse and transform data
end

stimulus = data[1][rat_num]             # Data frames obtained from a gaussian distribution of a determined coherence
coherences = data[2][rat_num]           # Target coherences used for each trial
rewards = data[3][rat_num]              # Correct choice side
decisions = data[4][rat_num]            # Actual decision made by the rat each trial
performances = data[5][rat_num]         # 0 if the rat chose the correct choice, 1 if the rat chose the incorrect choice
target_sigmas = data[6][rat_num]        # Target sigmas used for each trial
dates = data[7][rat_num]

stim_10 = []
stim_20 = []
choices_10 = []
choices_20 = []
coherences_10 = []
coherences_20 = []
rewards_10 = []
rewards_20 = []
sigma_10 = []
sigma_20 = []
performance_10 = []
performance_20 = []
date_10 = []
date_20 = []

# Separar els estímuls de 10 frames dels de 20 frames
for session in 1:size(stimulus)[1]
    stim_1 = stimulus[session]
    if size(stim_1)[1] == 10
        stim_10 = [stim_10;stim_1]
        choices_10 = [choices_10;decisions[session]]
        coherences_10 = [coherences_10;coherences[session]]
        rewards_10 = [rewards_10;rewards[session]]
        sigma_10 = [sigma_10;target_sigmas[session]]
        performance_10 = [performance_10;performances[session]]
        date_10 = [date_10;dates[session]]
    end
    if size(stim_1)[1] == 20
        stim_20 = [stim_20;stim_1]
        choices_20 = [choices_20;decisions[session]]
        coherences_20 = [coherences_20;coherences[session]]
        rewards_20 = [rewards_20;rewards[session]]
        sigma_20 = [sigma_20;target_sigmas[session]]
        performance_20 = [performance_20;performances[session]]
        date_20 = [date_20;dates[session]]
    end
end
n_stim_10 = convert(Int64, size(stim_10)[1]/10)
n_stim_20 = convert(Int64, size(stim_20)[1]/20)

results = Dict("stim_10"=>stim_10,"stim_20"=>stim_20,"choices_10"=>choices_10,"choices_20"=>choices_20,"n_stim_10"=>n_stim_10,"n_stim_20"=>n_stim_20,
"coherences_10"=>coherences_10,"coherences_20"=>coherences_20,"rewards_10"=>rewards_10,"rewards_20"=>rewards_20,"sigma_10"=>sigma_10,
"sigma_20"=>sigma_20,"performance_10"=>performance_10,"performance_20"=>performance_20,"date_10"=>date_10,"date_20"=>date_20)

filename_save = path_load*"data/processed_data_rat"*string(rat_num)*"_dataset2_1.json"
# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(results)

# write the file with the stringdata variable information
open(filename_save, "w") do f
        write(f, stringdata)
     end

max_stars_repo_name: emmafdz/FinalDegreeProject

max_issues_repo_name: emmafdz/FinalDegreeProject

max_forks_repo_name: emmafdz/FinalDegreeProject


text:
#!/usr/bin/env julia

using FourierFlows
using Test

# Run tests

time = @elapsed begin

println("-- Core tests --")

@testset "Grid tests" begin
  include("test_grid.jl")
end

@testset "FFT tests" begin
  include("test_fft.jl")
end

@testset "IFFT tests" begin
  include("test_ifft.jl")
end

@testset "Utils tests" begin
  include("test_utils.jl")
end

@testset "Timestepper tests" begin
  include("test_timesteppers.jl")
end

@testset "Diagnostics tests" begin
  include("test_diagnostics.jl")
end

println("-- Physics tests --")

@testset "Physics: Kuramoto-Sivashinsky" begin
  include("test_kuramotosivashinsky.jl")
end

@testset "Physics: TracerAdvDiff" begin
  include("test_traceradvdiff.jl")
end

end
println("Total test time: ", time)


text:
###############################################################################
#
#   Groups.jl : Nemo groups
#
###############################################################################

include("Rings.jl")

################################################################################
#
#  Element types for instances of groups
#
################################################################################

elem_type(::T) where {T <: Group} = elem_type(T)