text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

"""
	t(args::Array{Number,1}...)::Matrix

Return an *affine transformation Matrix* in homogeneous coordinates. Such `translation` Matrix has ``d+1`` rows and ``d+1`` columns, where ``d`` is the number of translation parameters in the `args` array.

# Examples

```julia
julia> t(1,2)			# 2D translation
# return
3×3 Array{Float64,2}:
 1.0  0.0  1.0
 0.0  1.0  2.0
 0.0  0.0  1.0

julia> Lar.t(1.,2,3)		# 3D translation
# return
4×4 Array{Float64,2}:
 1.0  0.0  0.0  1.0
 0.0  1.0  0.0  2.0
 0.0  0.0  1.0  3.0
 0.0  0.0  0.0  1.0
```
"""
function t(args...)
	d = length(args)
	mat = Matrix{Float64}(LinearAlgebra.I, d+1, d+1)
	for k in range(1, length=d)
        	mat[k,d+1]=args[k]
	end
	return mat
end

"""
	s(args::Array{Number,1}...)::Matrix

Return an *affine transformation Matrix* in homogeneous coordinates. Such `scaling` Matrix has ``d+1`` rows and ``d+1`` columns, where ``d`` is the number of scaling parameters in the `args` array.

# Examples

```julia
julia> Lar.s(2,3)			# 2D scaling
# return
3×3 Array{Float64,2}:
 2.0  0.0  0.0
 0.0  3.0  0.0
 0.0  0.0  1.0

julia> Lar.s(2.,3.,4.)		# 3D scaling
# return
4×4 Array{Float64,2}:
 2.0  0.0  0.0  0.0
 0.0  3.0  0.0  0.0
 0.0  0.0  4.0  0.0
 0.0  0.0  0.0  1.0
```
"""
function s(args...)
	d = length(args)
	mat = Matrix{Float64}(LinearAlgebra.I, d+1, d+1)
	for k in range(1, length=d)
		mat[k,k]=args[k]
	end
	return mat
end

"""
	r(args...)

Return an *affine transformation Matrix* in homogeneous coordinates. Such `Rotation` Matrix has *dimension* either equal to 3 or to 4, for 2D and 3D rotation, respectively.
The `{Number,1}` of `args` either contain a single `angle` parameter in *radiants*, or a vector with three elements, whose `norm` is the *rotation angle* in 3D and whose `normalized value` gives the direction of the *rotation axis* in 3D.

# Examples

```julia
julia> Lar.r(pi/6)				# 2D rotation of ``π/6`` angle
# return
3×3 Array{Float64,2}:
 0.866025  -0.5       0.0
 0.5        0.866025  0.0
 0.0        0.0       1.0

julia> Lar.r(0,0,pi/4)
# return
4×4 Array{Float64,2}:		# 3D rotation about the ``z`` axis, with ``π/6`` angle
 0.707107  -0.707107  0.0  0.0
 0.707107   0.707107  0.0  0.0
 0.0        0.0       1.0  0.0
 0.0        0.0       0.0  1.0

julia> Lar.r(1,1,1)		# 3D rotation about the ``x=y=z`` axis, with angle ``1.7320508`` angle
# return
4×4 Array{Float64,2}:
  0.226296  -0.183008   0.956712  0.0
  0.956712   0.226296  -0.183008  0.0
 -0.183008   0.956712   1.21332   0.0
  0.0        0.0        0.0       1.0

```
"""
function r(args...)
    args = collect(args)
    n = length(args)
    if n == 1 # rotation in 2D
        angle = args[1]; COS = cos(angle); SIN = sin(angle)
        mat = Matrix{Float64}(LinearAlgebra.I, 3, 3)
        mat[1,1] = COS;    mat[1,2] = -SIN;
        mat[2,1] = SIN;    mat[2,2] = COS;
    end

     if n == 3 # rotation in 3D
        mat = Matrix{Float64}(LinearAlgebra.I, 4, 4)
        angle = norm(args);
        if norm(args) != 0.0
			axis = args #normalize(args)
			COS = cos(angle); SIN= sin(angle)
			if axis[2]==axis[3]==0.0    # rotation about x
				mat[2,2] = COS;    mat[2,3] = -SIN;
				mat[3,2] = SIN;    mat[3,3] = COS;
			elseif axis[1]==axis[3]==0.0   # rotation about y
				mat[1,1] = COS;    mat[1,3] = SIN;
				mat[3,1] = -SIN;    mat[3,3] = COS;
			elseif axis[1]==axis[2]==0.0    # rotation about z
				mat[1,1] = COS;    mat[1,2] = -SIN;
				mat[2,1] = SIN;    mat[2,2] = COS;
			else
				I = Matrix{Float64}(LinearAlgebra.I, 3, 3); u = axis
				Ux=[0 -u[3] u[2] ; u[3] 0 -u[1] ;  -u[2] u[1] 1]
				UU =[u[1]*u[1]    u[1]*u[2]   u[1]*u[3];
					 u[2]*u[1]    u[2]*u[2]   u[2]*u[3];
					 u[3]*u[1]    u[3]*u[2]   u[3]*u[3]]
				mat[1:3,1:3]=COS*I+SIN*Ux+(1.0-COS)*UU
			end
		end
	end
	return mat
end

"""
	removeDups(CW::Cells)::Cells

Remove dublicate `cells` from `Cells` object. Then put `Cells` in *canonical form*, i.e. with *sorted indices* of vertices in each (unique) `Cells` Array element.
"""
function removeDups(CW::Cells)::Cells
	CW = collect(Set(CW))
	CWs = collect(map(sort,CW))
	return CWs
end

"""
	Struct

A *container* of geometrical objects is defined by applying the function `Struct` to
the array of contained objects. Each value is defined in local coordinates and may be transformed by affine transformation tensors.

The value returned from the application of `Struct` to an `Array{Any, 1}` of `LAR` values, `matrices`, and `Struct` values is a value of
`Struct type`.  The coordinate system of this value is the one associated with the first object of the `Struct` arguments.  Also,
the resulting geometrical value is often associated with a variable name.

The generation of containers may continue hierarchically by suitably applying `Struct`. Notice that each LAR object in a `Struct` container is transformed by each matrix before it *within the container*, going from right to left. The action of a transformation (tensor) extends to each object on its right within its own container. Whereas,  the action of a tensor does not extend outside its container, according to the semantics of *PHIGS* structures.

# Example

```julia
julia> L = LinearAlgebraicRepresentation;

julia> assembly = L.Struct([L.sphere()(), L.t(3,0,-1), L.cylinder()()])
# return
Lar.Struct(Any[([0.0 -0.173648 … -0.336824 -0.17101; 0.0 0.0 … 0.0593912 0.0301537;
-1.0 -0.984808 … 0.939693 0.984808], Array{Int64,1}[[2, 3, 1], [4, 2, 3], [4, 3, 5], [4,
5, 6], [7, 5, 6], [7, 8, 6], [7, 9, 8], … , [1.0 0.0 0.0 3.0; 0.0 1.0 0.0 0.0; 0.0 0.0 1.0
-1.0; 0.0 0.0 0.0 1.0], ([0.5 0.5 … 0.492404 0.492404; 0.0 0.0 … -0.0868241 -0.0868241;
0.0 2.0 … 0.0 2.0], Array{Int64,1}[[4, 2, 3, 1], [4, 3, 5, 6], [7, 5, 8, 6], [7, 9, 10,
8], [9, 10, 11, 12], [13, 14, 11, 12], … , [68, 66, 67, 65], [68, 69, 67, 70], [71, 69,
72, 70], [71, 2, 72, 1]])], Array{Float64,2}[[-1.0; -1.0; -1.0], [3.5; 1.0; 1.0]],
"14417445522513259426", 3, "feature")

julia> assembly.name = "simple example"
# return
"simple example"

julia> assembly
# return
Lar.Struct(Any[([0.0 -0.173648 … -0.336824 -0.17101; 0.0 0.0 … 0.0593912 0.0301537;
-1.0 -0.984808 … 0.939693 0.984808], Array{Int64,1}[[2, 3, 1], [4, 2, 3], [4, 3, 5], [4,
5, 6], [7, 5, 6], [7, 8, 6], … , [71, 2, 72, 1]])], Array{Float64,2}[[-1.0; -1.0; -1.0],
[3.5; 1.0; 1.0]], "simple example", 3, "feature")

julia> using Plasm

julia> Plasm.view(assembly)
```
"""
mutable struct Struct
	body::Array
	box
	name::AbstractString
	dim
	category::AbstractString

	function Struct()
		self = new([],Any,"new",Any,"feature")
		self.name = string(objectid(self))
		return self

	end

	function Struct(data::Array)
		self = Struct()
		self.body = data
		self.box = box(self)
		self.dim = length(self.box[1])
		return self
	end

	function Struct(data::Array,name)
		self = Struct()
		self.body=[item for item in data]
		self.box = box(self)
		self.dim = length(self.box[1])
		self.name = string(name)
		return self
	end

	function Struct(data::Array,name,category)
		self = Struct()
		self.body = [item for item in data]
		self.box = box(self)
		self.dim = length(self.box[1])
		self.name = string(name)
		self.category = string(category)
		return self
	end

end

	function name(self::Struct)
		return self.name
	end
	function category(self::Struct)
		return self.category
	end

	function len(self::Struct)
		return length(self.body)
	end
	function getitem(self::Struct,i::Int)
		return self.body[i]
	end
	function setitem(self::Struct,i,value)
		self.body[i]=value
	end
	function pprint(self::Struct)
		return "<Struct name: $(self.name)"
	end
	function set_name(self::Struct,name)
		self.name = string(name)
	end
	function clone(self::Struct,i=0)
		newObj = deepcopy(self)
		if i!=0
			newObj.name="$(self.name)_$(string(i))"
		end
		return newObj
	end
	function set_category(self::Struct,category)
		self.category = string(category)
	end

"""
	struct2lar(structure::Struct)::Union{LAR,LARmodel}

"""
function struct2lar(structure) # TODO: extend to true `LARmodels`
	listOfModels = Lar.evalStruct(structure)
	vertDict= Dict()
	index,defaultValue = 0,0
	W = Array{Float64,1}[]
	m = length(listOfModels[1])
	larmodel = [Array{Number,1}[] for k=1:m]

	for model in listOfModels
		V = model[1]
		for k=2:m
			for incell in model[k]
				outcell=[]
				for v in incell
					key = map(Lar.approxVal(7), V[:,v])
					if get(vertDict,key,defaultValue)==defaultValue
						index += 1
						vertDict[key]=index
						push!(outcell,index)
						push!(W,key)
					else
						push!(outcell,vertDict[key])
					end
				end
				append!(larmodel[k],[outcell])
			end
		end
	end

	append!(larmodel[1], W)
	V = hcat(larmodel[1]...)
	chains = [convert(Lar.Cells, chain) for chain in larmodel[2:end]]
	return (V, chains...)
end

"""
	embedTraversal(cloned::Struct,obj::Struct,n::Int,suffix::String)

# TODO:  debug embedTraversal
"""
function embedTraversal(cloned::Struct,obj::Struct,n::Int,suffix::String)
	for i=1:length(obj.body)
		if isa(obj.body[i],Matrix)
			mat = obj.body[i]
			d,d = size(mat)
			newMat = Matrix{Float64}(LinearAlgebra.I, d+n, d+n)
			for h in range(1, length=d)
				for k in range(1, length=d)
					newMat[h,k]=mat[h,k]
				end
			end
			push!(cloned.body,[newMat])
		elseif (isa(obj.body[i],Tuple) ||isa(obj.body[i],Array))&& length(obj.body[i])==3
			V,FV,EV = deepcopy(obj.body[i])
			dimadd = n
			ncols = size(V,2)
			nmat = zeros(dimadd,ncols)
			V = [V;zeros(dimadd,ncols)]
			push!(cloned.body,[(V,FV,EV)])
		elseif  (isa(obj.body[i],Tuple) ||isa(obj.body[i],Array))&& length(obj.body[i])==2
			V,EV = deepcopy(obj.body[i])
			dimadd = n
			ncols = size(V,2)
			nmat = zeros(dimadd,ncols)
			V = [V;zeros(dimadd,ncols)]
			push!(cloned.body,[(V,EV)])
		elseif isa(obj.body[i],Struct)
			newObj = Struct()
			newObj.box = [ [obj.body[i].box[1];zeros(dimadd)],
							[obj.body[i].box[2];zeros(dimadd)] ]
			newObj.category = (obj.body[i]).category
			push!(cloned.body,[embedTraversal(newObj,obj.body[i],obj.dim+n,suffix)])
		end
	end
	return cloned
end

"""
	embedStruct(n::Int)(self::Struct,suffix::String="New")

# TODO:  debug embedStruct
"""
function embedStruct(n::Int)
	function embedStruct0(self::Struct,suffix::String="New")
		if n==0
			return self, length(self.box[1])
		end
		cloned = Struct()
		cloned.box = [ [self.body[i].box[1];zeros(dimadd)],
						[self.body[i].box[2];zeros(dimadd)] ]
		cloned.name = self.name*suffix
		cloned.category = self.category
		cloned.dim = self.dim+n
		cloned = embedTraversal(cloned,self,n,suffix)
		return cloned
	end
	return embedStruct0
end

"""
	box(model)

"""
function box(model)
	if isa(model,Matrix)
		return nothing
	elseif isa(model,Struct)
		listOfModels = evalStruct(model)
		#dim = checkStruct(listOfModels)
		if listOfModels == []
			return model.box
		else
			theMin,theMax = box(listOfModels[1])
			for theModel in listOfModels[2:end]
				modelMin,modelMax= box(theModel)
				for (k,val) in enumerate(modelMin)
					if val < theMin[k]
						theMin[k]=val
					end
				end
				for (k,val) in enumerate(modelMax)
					if val > theMax[k]
						theMax[k]=val
					end
				end
			end
		end
		return [theMin,theMax]

	elseif (isa(model,Tuple) ||isa(model,Array))&& (length(model)>=2)
		V = model[1]
		theMin = minimum(V, dims=2)
		theMax = maximum(V, dims=2)
	end

	return [theMin,theMax]
end

"""
	apply(affineMatrix::Array{Float64,2}, larmodel::Union{LAR,LARmodel})
"""
function apply(affineMatrix, larmodel)
	data = collect(larmodel)
	V = data[1]

	m,n = size(V)
	W = [V; fill(1.0, (1,n))]
	V = (affineMatrix * W)[1:m,1:n]

	data[1] = V
	larmodel = Tuple(data)
	return larmodel
end

"""
	checkStruct(lst)

"""
function checkStruct(lst)
	obj = lst[1]
	if isa(obj,Matrix)
		dim = size(obj,1)-1
	elseif (isa(obj,Tuple) || isa(obj,Array))
		dim = length(obj[1][:,1])

	elseif isa(obj,Struct)
		dim = length(obj.box[1])
	end
	return dim
end

"""
	traversal(CTM,stack,obj,scene=[])

"""
function traversal(CTM::Matrix, stack, obj, scene=[])
	for i = 1:length(obj.body)
		if isa(obj.body[i],Matrix)
			CTM = CTM*obj.body[i]
		elseif (isa(obj.body[i],Tuple) || isa(obj.body[i],Array)) &&
			(length(obj.body[i])>=2)
			l = apply(CTM, obj.body[i])
			push!(scene,l)
		elseif isa(obj.body[i],Struct)
			push!(stack,CTM)
			traversal(CTM,stack,obj.body[i],scene)
			CTM = pop!(stack)
		end
	end
	return scene
end

"""
	evalStruct(self)

"""
function evalStruct(self::Struct)
	dim = checkStruct(self.body)
   	CTM, stack = Matrix{Float64}(LinearAlgebra.I, dim+1, dim+1), []
   	scene = traversal(CTM, stack, self, [])
	return scene
end


text:
using Flux, Test, Mill
using Mill: reflectinmodel
using Combinatorics

# use only activations without "kinks" for numerical gradient checking
# see e.g. https://stackoverflow.com/questions/40623512/how-to-check-relu-gradient
const ACTIVATIONS = [identity, σ, swish, softplus, logcosh, mish, tanhshrink, lisht]

@testset "aggregation grad check w.r.t. input" begin
    for bags in BAGS2
        d = rand(1:20)
        x = randn(d, 10)
        w = abs.(randn(size(x, 2))) .+ 0.1
        w_mat = abs.(randn(size(x))) .+ 0.1

        # generate all combinations of aggregations
        anames = ["Sum", "Mean", "Max", "PNorm", "LSE"]
        for idxs in powerset(collect(1:length(anames)))
            !isempty(idxs) || continue
            # not a thorough testing of all functions, but fast enough
            length(idxs) <= 3 || continue
            # for idxs in permutations(idxs)

            s = Symbol("Segmented", anames[idxs]...)
            a = @eval $s($d)
            @test mgradtest(x) do x
                a(x, bags)
            end
            @test mgradtest(x) do x
                a(x, bags, w)
            end
            @test mgradtest(x) do x
                a(x, bags, w_mat)
            end
        end
    end
end

@testset "aggregation grad check w.r.t. agg params" begin
    # r_map and p_map are stable
    @test first(gradient(softplus, 10000)) ≈ σ(10000) ≈ 1.0
    @test first(gradient(softplus, -10000)) ≈ σ(-10000) ≈ 0

    fs = [:SegmentedSum, :SegmentedMax, :SegmentedMean, :SegmentedPNorm, :SegmentedLSE]
    params = [(:C1,), (:C2,), (:C3,), (:ρ1, :c, :C4), (:ρ2, :C5)]

    for idxs in powerset(collect(1:length(fs)))
        !isempty(idxs) || continue;
        length(idxs) <= 2 || continue

        d = rand(1:20)
        x = randn(d, 0)
        as = []; cs = []; rs = []
        for (f, ps) in zip(fs[idxs], params[idxs])
            push!(rs, fill(:(randn($d)), length(ps))...)
            push!(as, ps...)
            push!(cs, Expr(:call, f, ps...))
        end
        @eval begin
            @test mgradtest($(map(eval, rs)...)) do $(as...)
                a = Aggregation($(cs...))
                a(missing, ScatteredBags([Int[], Int[]]))
            end
            @test mgradtest($(map(eval, rs)...)) do $(as...)
                a = Aggregation($(cs...))
                a(missing, AlignedBags([0:-1]), nothing)
            end
            @test mgradtest($(map(eval, rs)...)) do $(as...)
                a = Aggregation($(cs...))
                a($x, ScatteredBags([Int[]]))
            end
            @test mgradtest($(map(eval, rs)...)) do $(as...)
                a = Aggregation($(cs...))
                a($x, AlignedBags([0:-1, 0:-1]), nothing)
            end
        end

        for bags in BAGS2
            d = rand(1:20)
            x = randn(d, 10)
            w = abs.(randn(size(x, 2))) .+ 0.1
            w_mat = abs.(randn(size(x))) .+ 0.1
            as = []; cs = []; rs = []
            for (f, ps) in zip(fs[idxs], params[idxs])
                push!(rs, fill(:(randn($d)), length(ps))...)
                push!(as, ps...)
                push!(cs, Expr(:call, f, ps...))
            end
            @eval begin
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a($x, $bags)
                end
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a($x, $bags, $w)
                end
                @test mgradtest($(map(eval, rs)...)) do $(as...)
                    a = Aggregation($(cs...))
                    a($x, $bags, $w_mat)
                end
            end
        end
    end
end

@testset "derivative w.r.t weights in aggregations" begin
    for bags in BAGS2
        d = rand(1:20)
        x = randn(d, 10)
        w = abs.(randn(size(x, 2))) .+ 0.1
        w_mat = abs.(randn(size(x))) .+ 0.1

        a1 = f64(SegmentedSum(d))
        a2 = f64(SegmentedMean(d))
        a3 = f64(SegmentedMax(d))
        a4 = f64(SegmentedPNorm(d))
        a5 = f64(SegmentedLSE(d))
        for g in [
                  w -> a1(x, bags, w),
                  w -> a2(x, bags, w),
                  w -> a3(x, bags, w),
                  w -> a5(x, bags, w)
                 ]
            @test mgradtest(g, w)
        end
        for g in [
                  w_mat -> a1(x, bags, w_mat),
                  w_mat -> a2(x, bags, w_mat),
                  w_mat -> a3(x, bags, w_mat),
                  w_mat -> a5(x, bags, w_mat)
                 ]
            @test mgradtest(g, w_mat)
        end
        # for g in [
        #           w -> a4(x, bags, w_mat)
        #          ]
        #     # NOT IMPLEMENTED YET
        #     @test_throws Exception mgradtest(g, w_mat)
        # end
        # for g in [
        #           w -> a4(x, bags, w_mat)
        #          ]
        #     # NOT IMPLEMENTED YET
        #     @test_throws Exception mgradtest(g, w_mat)
        # end
    end
end

@testset "model aggregation grad check w.r.t. inputs" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2, rand(ACTIVATIONS))
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)

        n = ArrayNode(x)
        m = f64(reflectinmodel(n, layerbuilder))
        @test mgradtest(x) do x
            n = ArrayNode(x)
            m(n).data
        end

        bn = BagNode(ArrayNode(x), bags1)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder))
        @test mgradtest(x) do x
            bn = BagNode(ArrayNode(x), bags1)
            m(bn).data
        end

        tn = ProductNode((ArrayNode(x), ArrayNode(y)))
        m = f64(reflectinmodel(tn, layerbuilder))
        @test mgradtest(x, y) do x, y
            tn = ProductNode((ArrayNode(x), ArrayNode(y)))
            m(tn).data
        end

        tn = ProductNode((BagNode(ArrayNode(y), bags1), BagNode(ArrayNode(x), bags2)))
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        @test mgradtest(x, y) do x, y
            tn = ProductNode((BagNode(ArrayNode(y), bags1), BagNode(ArrayNode(x), bags2)))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3)
        bnn = BagNode(bn, bags1)
        abuilder = d -> SegmentedPNormLSESumMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        @test mgradtest(z) do z
            bn = BagNode(ArrayNode(z), bags3)
            bnn = BagNode(bn, bags1)
            m(bnn).data
        end
    end
end

@testset "model aggregation grad check w.r.t. inputs weighted" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2, rand(ACTIVATIONS))
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)
        w = abs.(randn(4)) .+ 0.1
        w2 = abs.(randn(4)) .+ 0.1
        w3 = abs.(randn(8)) .+ 0.1

        bn = BagNode(ArrayNode(x), bags1, w)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder))
        @test mgradtest(x) do x
            bn = BagNode(ArrayNode(x), bags1, w)
            m(bn).data
        end

        tn = ProductNode((BagNode(ArrayNode(y), bags1, w), BagNode(ArrayNode(x), bags2, w2)))
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        @test mgradtest(x, y) do x, y
            tn = ProductNode((BagNode(ArrayNode(y), bags1, w), BagNode(ArrayNode(x), bags2, w2)))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3, w3)
        bnn = BagNode(bn, bags1)
        abuilder = d -> SegmentedPNormLSESumMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        @test mgradtest(z) do z
            bn = BagNode(ArrayNode(z), bags3, w3)
            bnn = BagNode(bn, bags1, w)
            m(bnn).data
        end
    end
end

@testset "model aggregation grad check w.r.t. params" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2)
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)

        n = ArrayNode(x)
        m = f64(reflectinmodel(n, layerbuilder))
        a = rand(ACTIVATIONS)
        @test mgradtest(params(m)...) do W, b
            m = ArrayModel(Dense(W, b, a))
            m(n).data
        end

        bn = BagNode(ArrayNode(x), bags1)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder, abuilder))
        a1, a2 = rand(ACTIVATIONS, 2)
        @test mgradtest(params(m)...) do W1, b1, ρ1, c, C1, ρ2, C2, W2, b2
            m = BagModel(Dense(W1, b1, a1),
                         Aggregation(
                                     SegmentedPNorm(ρ1, c, C1),
                                     SegmentedLSE(ρ2, C2)
                                    ),
                         Dense(W2, b2, a2))
            m(bn).data
        end

        tn = ProductNode((ArrayNode(x), ArrayNode(y)))
        m = f64(reflectinmodel(tn, layerbuilder))
        a1, a2, a3 = rand(ACTIVATIONS, 3)
        @test mgradtest(params(m)...) do W1, b1, W2, b2, W3, b3
            m = ProductModel(ArrayModel.((
                                          Dense(W1, b1, a1),
                                          Dense(W2, b2, a2)
                                         )), Dense(W3, b3, a3))
            m(tn).data
        end

        tn = ProductNode((BagNode(ArrayNode(y), bags1), BagNode(ArrayNode(x), bags2)))
        abuilder = d -> SegmentedPNormLSESumMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        a1, a2, a3, a4, a5 = rand(ACTIVATIONS, 5)
        @test mgradtest(params(m)...) do W1, b1, ρ11, c1, C11, ρ12, C12, C13, C14,
            W2, b2, W3, b3, ρ21, c2, C21, ρ22, C22, C23, C24, W4, b4, W5, b5
            m = ProductModel((
                              BagModel(
                                       Dense(W1, b1, a1),
                                       Aggregation(
                                                   SegmentedPNorm(ρ11, c1, C11),
                                                   SegmentedLSE(ρ12, C12),
                                                   SegmentedSum(C13),
                                                   SegmentedMax(C14)
                                                  ),
                                       Dense(W2, b2, a2)
                                      ),
                              BagModel(
                                       Dense(W3, b3, a3),
                                       Aggregation(
                                                   SegmentedPNorm(ρ21, c2, C21),
                                                   SegmentedLSE(ρ22, C22),
                                                   SegmentedSum(C23),
                                                   SegmentedMax(C24)
                                                  ),
                                       Dense(W4, b4, a4)
                                      ),
                             ), Dense(W5, b5, a5))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3)
        bnn = BagNode(bn, bags1)
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        a1, a2, a3 = rand(ACTIVATIONS, 3)
        @test mgradtest(params(m)...) do W1, b1, C11, C12, W2, b2, C21, C22, W3, b3
            m = BagModel(
                         BagModel(
                                  Dense(W1, b1, a1),
                                  Aggregation(
                                              SegmentedMean(C11),
                                              SegmentedMax(C12)
                                             ),
                                  Dense(W2, b2, a2)
                                 ),
                         Aggregation(
                                     SegmentedMean(C21),
                                     SegmentedMax(C22)
                                    ),
                         Dense(W3, b3, a3)
                        )
            m(bnn).data
        end
    end
end

@testset "model aggregation grad check w.r.t. params weighted" begin
    for (bags1, bags2, bags3) in BAGS3
        layerbuilder(k) = Dense(k, 2)
        x = randn(4, 4)
        y = randn(3, 4)
        z = randn(2, 8)
        w = abs.(randn(4)) .+ 0.1
        w2 = abs.(randn(4)) .+ 0.1
        w3 = abs.(randn(8)) .+ 0.1

        bn = BagNode(ArrayNode(x), bags1, w)
        abuilder = d -> SegmentedPNormLSE(d)
        m = f64(reflectinmodel(bn, layerbuilder, abuilder))
        a1, a2 = rand(ACTIVATIONS, 2)
        @test mgradtest(params(m)...) do W1, b1, ρ1, c, C1, ρ2, C2, W2, b2
            m = BagModel(Dense(W1, b1, a1),
                         Aggregation(
                                     SegmentedPNorm(ρ1, c, C1),
                                     SegmentedLSE(ρ2, C2)
                                    ),
                         Dense(W2, b2, a2))
            m(bn).data
        end

        tn = ProductNode((BagNode(ArrayNode(y), bags1, w), BagNode(ArrayNode(x), bags2, w2)))
        abuilder = d -> SegmentedPNormLSESumMax(d)
        m = f64(reflectinmodel(tn, layerbuilder, abuilder))
        a1, a2, a3, a4, a5 = rand(ACTIVATIONS, 5)
        @test mgradtest(params(m)...) do W1, b1, ρ11, c1, C11, ρ12, C12, C13, C14,
            W2, b2, W3, b3, ρ21, c2, C21, ρ22, C22, C23, C24, W4, b4, W5, b5
            m = ProductModel((
                              BagModel(
                                       Dense(W1, b1, a1),
                                       Aggregation(
                                                   SegmentedPNorm(ρ11, c1, C11),
                                                   SegmentedLSE(ρ12, C12),
                                                   SegmentedSum(C13),
                                                   SegmentedMax(C14)
                                                  ),
                                       Dense(W2, b2, a2)
                                      ),
                              BagModel(
                                       Dense(W3, b3, a3),
                                       Aggregation(
                                                   SegmentedPNorm(ρ21, c2, C21),
                                                   SegmentedLSE(ρ22, C22),
                                                   SegmentedSum(C23),
                                                   SegmentedMax(C24)
                                                  ),
                                       Dense(W4, b4, a4)
                                      ),
                             ), Dense(W5, b5, a5))
            m(tn).data
        end

        bn = BagNode(ArrayNode(z), bags3, w3)
        bnn = BagNode(bn, bags1, w)
        abuilder = d -> SegmentedMeanMax(d)
        m = f64(reflectinmodel(bnn, layerbuilder, abuilder))
        a1, a2, a3 = rand(ACTIVATIONS, 3)
        @test mgradtest(params(m)...) do W1, b1, C11, C12, W2, b2, C21, C22, W3, b3
            m = BagModel(
                         BagModel(
                                  Dense(W1, b1, a1),
                                  Aggregation(
                                              SegmentedMean(C11),
                                              SegmentedMax(C12)
                                             ),
                                  Dense(W2, b2, a2)
                                 ),
                         Aggregation(
                                     SegmentedMean(C21),
                                     SegmentedMax(C22)
                                    ),
                         Dense(W3, b3, a3)
                        )
            m(bnn).data
        end
    end

    @testset "A gradient of ProductNode with a NamedTuple " begin
        a1, a2, a3 = rand(ACTIVATIONS, 3)
        x = ProductNode((
                         a = ArrayNode(randn(2, 4)),
        b = ArrayNode(randn(3, 4))
       ))
        m = f64(ProductModel((
                              a = ArrayModel(Dense(2, 2, a1)),
                              b = ArrayModel(Dense(3, 1, a2))
                             ), Dense(3, 2, a3)))
        @test mgradtest(params(m)...) do W1, b1, W2, b2, W3, b3
            m = ProductModel((
                              a = ArrayModel(Dense(W1, b1, a1)),
            b = ArrayModel(Dense(W2, b2, a2))
           ), Dense(W3, b3, a3))
            sum(m(x).data)
        end
    end
end

ext: jl


text:
__precompile__(true)

"""
The PeriodicTable module exports an `elements` variable that stores
data (of type `Element`) on every element in the periodic table.

The data can be looked up by atomic number via `elements[n]`, by name
(case-insensitive) via e.g. `elements["oxygen"]`, or by symbol via
e.g. `elements[:O]`.
"""
module PeriodicTable
export Element, elements

import Unitful: u, g, cm, K, J, mol, Quantity

"""
Element composite type
"""
struct Element
    name::String
    appearance::String
    atomic_mass::typeof(1.0u)
    boil::typeof(1.0K)
    category::String
    color::String
    density::typeof(1.0g/cm^3)
    discovered_by::String
    el_config::String
    melt::typeof(1.0K)
    molar_heat::typeof(1.0J/(mol*K))
    named_by::String
    number::Int
    period::Int
    phase::String
    source::String # url
    spectral_img::String # url
    summary::String
    symbol::String
    xpos::Int
    ypos::Int
    shells::Vector{Int}
end

Element(; name::AbstractString="",
          appearance::AbstractString="",
          atomic_mass::typeof(1.0u)=NaN*u,
          boil::typeof(1.0K)=NaN*K,
          category::AbstractString="",
          color::AbstractString="",
          density::typeof(1.0g/cm^3)=NaN*g/cm^3,
          discovered_by::AbstractString="",
          el_config::AbstractString="",
          melt::typeof(1.0K)=NaN*K,
          molar_heat::typeof(1.0J/(mol*K))=NaN*J/(mol*K),
          named_by::AbstractString="",
          number::Integer=-1,
          period::Integer=-1,
          phase::AbstractString="",
          source::AbstractString="",
          spectral_img::AbstractString="",
          summary::AbstractString="",
          symbol::AbstractString="",
          xpos::Integer=-1,
          ypos::Integer=-1,
          shells::AbstractVector{<:Integer}=Int[]) =
    Element(name, appearance, atomic_mass, boil, category, color, density,
        discovered_by, el_config, melt, molar_heat, named_by, number, period, phase,
        source, spectral_img, summary, symbol, xpos, ypos, shells)

Base.show(io::IO, el::Element) = print(io, "Element(", el.name, ')')

ispresent(s) = !isempty(s)
ispresent(x::Union{Float64, Quantity}) = !isnan(x)
ispresent(n::Int) = n ≥ 0
function printpresent(io::IO, name, val, suffix=""; pad=16)
    if ispresent(val)
        println(io, lpad(name, pad), ": ", typeof(val) <: Quantity ? val.val : val, suffix)
    end
end

function Base.show(io::IO, ::MIME"text/plain", el::Element)
    println(io, el.name, " (", el.symbol, "), number ", el.number, ':')
    printpresent(io, "category", el.category)
    printpresent(io, "atomic mass", el.atomic_mass, " u")
    printpresent(io, "density", el.density, " g/cm³")
    printpresent(io, "molar heat", el.molar_heat, " J/mol⋅K")
    printpresent(io, "melting point", el.melt, " K")
    printpresent(io, "boiling point", el.boil, " K")
    printpresent(io, "phase", el.phase)
    printpresent(io, "shells", el.shells)
    printpresent(io, "e⁻-configuration", el.el_config)
    printpresent(io, "appearance", el.appearance)
    printpresent(io, "color", el.color)
    printpresent(io, "summary", el.summary)
    printpresent(io, "discovered by", el.discovered_by)
    printpresent(io, "named by", el.named_by)
    printpresent(io, "source", el.source)
    printpresent(io, "spectral image", el.spectral_img)
end

function printpresenthtml(io::IO, name, val, suffix="")
    if ispresent(val)
        println(io, "
", name, "	", typeof(val) <: Quantity ? val.val : val, suffix, "
")
    end
end

function Base.show(io::IO, ::MIME"text/html", el::Element)
    println(io, "")
    println(io, el.name, " (", el.symbol, "), number ", el.number, ':')
    println(io, "
")
    printpresenthtml(io, "category", el.category)
    printpresenthtml(io, "atomic mass", el.atomic_mass, " u")
    printpresenthtml(io, "density", el.density, " g/cm³")
    printpresenthtml(io, "molar heat", el.molar_heat, " J/mol⋅K")
    printpresenthtml(io, "melting point", el.melt, " K")
    printpresenthtml(io, "boiling point", el.boil, " K")
    printpresenthtml(io, "phase", el.phase)
    printpresenthtml(io, "shells", el.shells)
    printpresenthtml(io, "electron configuration", el.el_config)
    printpresenthtml(io, "appearance", el.appearance)
    printpresenthtml(io, "color", el.color)
    printpresenthtml(io, "summary", el.summary)
    printpresenthtml(io, "discovered by", el.discovered_by)
    printpresenthtml(io, "named by", el.named_by)

    link = string("", el.source, "")
    printpresenthtml(io, "source", link)
    println(io, "
")

    if ispresent(el.spectral_img)
        width = 500
        file = replace(el.spectral_img, "https://en.wikipedia.org/wiki/File:" => "")
        #Wikimedia hotlinking api
        wm = "https://commons.wikimedia.org/w/index.php?title=Special:Redirect/file/"
        imgdomain = string(wm, file, "&width=$width")
        println(io, "")
    end
end

"""
Elements composite type
"""
struct Elements
    data::Vector{Element}
    bynumber::Dict{Int,Element}
    byname::Dict{String,Element}
    bysymbol::Dict{Symbol,Element}
    Elements(data::Vector{Element}) = new(
        sort!(data, by=d->d.number),
        Dict{Int,Element}(d.number=>d for d in data),
        Dict{String,Element}(lowercase(d.name)=>d for d in data),
        Dict{Symbol,Element}(Symbol(d.symbol)=>d for d in data))
end
Base.getindex(e::Elements, i::Integer) = e.bynumber[i]
Base.getindex(e::Elements, i::AbstractString) = e.byname[lowercase(i)]
Base.getindex(e::Elements, i::Symbol) = e.bysymbol[i]
Base.getindex(e::Elements, v::AbstractVector) = Element[e[i] for i in v]
Base.haskey(e::Elements, i::Integer) = haskey(e.bynumber, i)
Base.haskey(e::Elements, i::AbstractString) = haskey(e.byname, lowercase(i))
Base.haskey(e::Elements, i::Symbol) = haskey(e.bysymbol, i)
Base.get(e::Elements, i::Integer, default) = get(e.bynumber, i, default)
Base.get(e::Elements, i::AbstractString, default) = get(e.byname, lowercase(i), default)
Base.get(e::Elements, i::Symbol, default) = get(e.bysymbol, i, default)

# support iterating over elements
Base.eltype(e::Elements) = Element
Base.length(e::Elements) = length(e.data)
if VERSION < v"0.7-"
    Base.start(e::Elements) = start(e.data)
    Base.next(e::Elements, i) = next(e.data, i)
    Base.done(e::Elements, i) = done(e.data, i)
else
    Base.iterate(e::Elements, state...) = iterate(e.data, state...)
end

# compact one-line printing
Base.show(io::IO, e::Elements) = print(io, "Elements(…", length(e), " elements…)")

# pretty-printing as a periodic table
function Base.show(io::IO, ::MIME"text/plain", e::Elements)
    println(io, e, ':')
    table = fill("   ", 10,18)
    for el in e
        table[el.ypos, el.xpos] = rpad(el.symbol, 3)
    end
    for i = 1:size(table,1)
        for j = 1:size(table, 2)
            print(io, table[i,j])
        end
        println(io)
    end
end

# Since Element equality is determined by atomic number alone...
Base.isequal(elm1::Element, elm2::Element) = elm1.number == elm2.number

# There is no need to use all the data in Element to calculated the hash
# since Element equality is determined by atomic number alone.
Base.hash(elm::Element, h::UInt) = hash(elm.number, h)

# Compare elements by atomic number to produce the most common way elements
# are sorted.
Base.isless(elm1::Element, elm2::Element) = elm1.number < elm2.number

# Provide a simple way to iterate over all elements.
Base.eachindex(elms::Elements) = eachindex(elms.data)

include("elements.jl")
const elements = Elements(_elements_data)

end


text:
using BlockArrays
using LinearAlgebra
using PredictiveControl

# Test a matrix kron into a matrix
A = ones( Float64, (2, 2) )
B = rand( 2, 2 )

K = blockkron( A, B )

@test blocksize( K ) == ( 2, 2 )
@test K[Block(1, 1)] == B
@test K[Block(1, 2)] == B
@test K[Block(2, 1)] == B
@test K[Block(2, 2)] == B

# Test a matrix kron into a vector
A = ones( Float64, (2) )
B = rand( 2, 2 )

K = blockkron( A, B )

@test blocksize( K ) == ( 2, 1 )
@test K[Block(1, 1)] == B
@test K[Block(2, 1)] == B

# Test a vector kron into a matrix
A = ones( Float64, (2, 2) )
B = rand(4,1)

K = blockkron( A, B )

@test blocksize( K ) == ( 2, 2 )
@test K[Block(1, 1)] == B
@test K[Block(1, 2)] == B
@test K[Block(2, 1)] == B
@test K[Block(2, 2)] == B

# Test a vector kron into a vector
A = ones( Float64, (2) )
B = rand(4,1)

K = blockkron( A, B )

@test blocksize( K ) == ( 2, 1 )
@test K[Block(1, 1)] == B
@test K[Block(2, 1)] == B

# Test a vector kron into a vector
A = ones( Float64, (2) )
B = rand(4)

K = blockkron( A, B )

@test blocksize( K ) == ( 2, )
@test reshape( K[Block(1, 1)], (:) ) == B
@test reshape( K[Block(2, 1)], (:) ) == B


text:
include("TwoPartCost.jl")
include("ThreePartCost.jl")
include("MarketBidCost.jl")
include("MultiStartCost.jl")
include("Area.jl")
include("LoadZone.jl")
include("Bus.jl")
include("Arc.jl")
include("Line.jl")
include("MonitoredLine.jl")
include("PhaseShiftingTransformer.jl")
include("TapTransformer.jl")
include("Transformer2W.jl")
include("HVDCLine.jl")
include("VSCDCLine.jl")
include("InterruptibleLoad.jl")
include("FixedAdmittance.jl")
include("PowerLoad.jl")
include("HydroEnergyReservoir.jl")
include("HydroDispatch.jl")
include("HydroPumpedStorage.jl")
include("RenewableDispatch.jl")
include("RenewableFix.jl")
include("ThermalStandard.jl")
include("ThermalMultiStart.jl")
include("GenericBattery.jl")
include("BatteryEMS.jl")
include("StaticReserve.jl")
include("StaticReserveNonSpinning.jl")
include("StaticReserveGroup.jl")
include("ReserveDemandCurve.jl")
include("VariableReserve.jl")
include("VariableReserveNonSpinning.jl")
include("AGC.jl")
include("Transfer.jl")
include("AVRFixed.jl")
include("AVRSimple.jl")
include("ESDC1A.jl")
include("ESDC2A.jl")
include("IEEET1.jl")
include("AVRTypeI.jl")
include("AVRTypeII.jl")
include("SCRX.jl")
include("ESAC1A.jl")
include("EXAC1A.jl")
include("EXAC1.jl")
include("EXAC2.jl")
include("ESAC6A.jl")
include("ESST1A.jl")
include("EXPIC1.jl")
include("ESST4B.jl")
include("BaseMachine.jl")
include("RoundRotorMachine.jl")
include("SalientPoleMachine.jl")
include("AndersonFouadMachine.jl")
include("FullMachine.jl")
include("MarconatoMachine.jl")
include("OneDOneQMachine.jl")
include("SimpleAFMachine.jl")
include("SimpleFullMachine.jl")
include("SimpleMarconatoMachine.jl")
include("PSSFixed.jl")
include("PSSSimple.jl")
include("IEEEST.jl")
include("SingleMass.jl")
include("FiveMassShaft.jl")
include("TGFixed.jl")
include("GasTG.jl")
include("GeneralGovModel.jl")
include("SteamTurbineGov1.jl")
include("HydroTurbineGov.jl")
include("IEEETurbineGov1.jl")
include("TGTypeI.jl")
include("TGTypeII.jl")
include("AverageConverter.jl")
include("FixedDCSource.jl")
include("ZeroOrderBESS.jl")
include("LCLFilter.jl")
include("LCFilter.jl")
include("KauraPLL.jl")
include("VirtualInertia.jl")
include("ReactivePowerDroop.jl")
include("CurrentControl.jl")
include("Source.jl")

export get_A1
export get_A2
export get_A3
export get_A4
export get_A5
export get_A6
export get_AT
export get_A_set
export get_Ae
export get_At
export get_Be
export get_D
export get_DB_h
export get_DB_l
export get_D_12
export get_D_23
export get_D_34
export get_D_45
export get_D_T
export get_D_ex
export get_D_hp
export get_D_ip
export get_D_lp
export get_D_turb
export get_Dm
export get_E_sat
export get_Efd_lim
export get_H
export get_H_ex
export get_H_hp
export get_H_ip
export get_H_lp
export get_I_lr
export get_K
export get_K0
export get_K1
export get_K2
export get_K3
export get_K4
export get_K5
export get_K6
export get_K7
export get_K8
export get_K_d
export get_K_ex
export get_K_hp
export get_K_i
export get_K_im
export get_K_ip
export get_K_ir
export get_K_lp
export get_K_lr
export get_K_p
export get_K_pm
export get_K_pr
export get_K_turb
export get_K_ω
export get_Ka
export get_Kb
export get_Kc
export get_Kd
export get_Kd_gov
export get_Ke
export get_Kf
export get_Kg
export get_Kh
export get_Ki
export get_Ki_gov
export get_Ki_load
export get_Ki_mw
export get_Kl
export get_Kp
export get_Kp_gov
export get_Kp_load
export get_Ks
export get_Kt
export get_Kv
export get_L_1d
export get_L_1q
export get_L_ad
export get_L_aq
export get_L_d
export get_L_f1d
export get_L_ff
export get_L_q
export get_Ld_ref
export get_Ls_lim
export get_PSS_flags
export get_P_ref
export get_R
export get_R_1d
export get_R_1q
export get_R_close
export get_R_f
export get_R_lim
export get_R_open
export get_R_th
export get_Rselect
export get_Se
export get_T1
export get_T2
export get_T3
export get_T4
export get_T5
export get_T6
export get_T7
export get_T_AA
export get_T_act
export get_T_eng
export get_T_rate
export get_Ta
export get_Ta_2
export get_Ta_3
export get_Ta_4
export get_Ta_Tb
export get_Tb
export get_Tb1
export get_Tc
export get_Tc1
export get_Td0_p
export get_Td0_pp
export get_Td_gov
export get_Te
export get_Tf
export get_Tf_1
export get_Tf_2
export get_Tf_load
export get_Tg
export get_Th
export get_Tj
export get_Tk
export get_Tpelec
export get_Tq0_p
export get_Tq0_pp
export get_Tr
export get_Ts
export get_Tsa
export get_Tsb
export get_Tw
export get_U0
export get_UEL_flags
export get_U_c
export get_VB_max
export get_VELM
export get_VFE_lim
export get_VH_max
export get_V_lim
export get_V_lr
export get_V_pss
export get_V_ref
export get_Va_lim
export get_Vcl
export get_Vcu
export get_Vdc_ref
export get_Vf
export get_Vi_lim
export get_Vm_lim
export get_Vr_lim
export get_Wf_nl
export get_X_th
export get_Xd
export get_Xd_p
export get_Xd_pp
export get_Xl
export get_Xq
export get_Xq_p
export get_Xq_pp
export get_Y
export get_active_power
export get_active_power_flow
export get_active_power_limits
export get_active_power_limits_from
export get_active_power_limits_pump
export get_active_power_limits_to
export get_ancillary_services
export get_angle
export get_angle_limits
export get_arc
export get_area
export get_available
export get_b
export get_base_power
export get_base_voltage
export get_battery_resistance
export get_battery_voltage
export get_bias
export get_bus
export get_bustype
export get_cf
export get_contributing_services
export get_conversion_factor
export get_db
export get_dc_dc_inductor
export get_dc_link_capacitance
export get_delta_t
export get_dynamic_injector
export get_efficiency
export get_energy_value
export get_eq_p
export get_ext
export get_fixed
export get_flow_limits
export get_from
export get_fs
export get_fuel
export get_fuel_flag
export get_gate_position_limits
export get_inflow
export get_initial_ace
export get_initial_energy
export get_initial_storage
export get_input_active_power_limits
export get_input_code
export get_internal_angle
export get_internal_voltage
export get_inv_d_fluxlink
export get_inv_q_fluxlink
export get_inverter_firing_angle
export get_inverter_tap_limits
export get_inverter_xrc
export get_kad
export get_kd
export get_kffi
export get_kffv
export get_ki_pll
export get_kic
export get_kii
export get_kiv
export get_kp_pll
export get_kpc
export get_kpi
export get_kpv
export get_kq
export get_kω
export get_lf
export get_lg
export get_load_response
export get_load_zone
export get_loss
export get_lv
export get_magnitude
export get_max_active_power
export get_max_reactive_power
export get_model
export get_must_run
export get_n_states
export get_no_load
export get_number
export get_operation_cost
export get_outflow
export get_output_active_power_limits
export get_peak_active_power
export get_peak_reactive_power
export get_penalty_cost
export get_power_factor
export get_power_trajectory
export get_primary_shunt
export get_prime_mover
export get_pump_efficiency
export get_q_nl
export get_r
export get_ramp_limits
export get_ramp_limits_pump
export get_rate
export get_rated_current
export get_rated_voltage
export get_rating
export get_rating_pump
export get_rc_rfd
export get_reactive_power
export get_reactive_power_flow
export get_reactive_power_limits
export get_reactive_power_limits_from
export get_reactive_power_limits_pump
export get_reactive_power_limits_to
export get_rectifier_firing_angle
export get_rectifier_tap_limits
export get_rectifier_xrc
export get_remote_bus_control
export get_requirement
export get_rf
export get_rg
export get_rv
export get_saturation_coeffs
export get_services
export get_shut_down
export get_speed_error_signal
export get_start_time_limits
export get_start_types
export get_start_up
export get_state_of_charge_limits
export get_states
export get_states_types
export get_status
export get_storage_capacity
export get_storage_target
export get_switch
export get_tap
export get_time_at_status
export get_time_frame
export get_time_limits
export get_time_limits_pump
export get_to
export get_valve_position_limits
export get_variable
export get_voltage
export get_voltage_limits
export get_x
export get_α
export get_γ_d1
export get_γ_d2
export get_γ_q1
export get_γ_q2
export get_γ_qd
export get_γd
export get_γq
export get_θp
export get_θp_rad
export get_τ_limits
export get_ω_lp
export get_ωad
export get_ωf
export set_A1!
export set_A2!
export set_A3!
export set_A4!
export set_A5!
export set_A6!
export set_AT!
export set_A_set!
export set_Ae!
export set_At!
export set_Be!
export set_D!
export set_DB_h!
export set_DB_l!
export set_D_12!
export set_D_23!
export set_D_34!
export set_D_45!
export set_D_T!
export set_D_ex!
export set_D_hp!
export set_D_ip!
export set_D_lp!
export set_D_turb!
export set_Dm!
export set_E_sat!
export set_Efd_lim!
export set_H!
export set_H_ex!
export set_H_hp!
export set_H_ip!
export set_H_lp!
export set_I_lr!
export set_K!
export set_K0!
export set_K1!
export set_K2!
export set_K3!
export set_K4!
export set_K5!
export set_K6!
export set_K7!
export set_K8!
export set_K_d!
export set_K_ex!
export set_K_hp!
export set_K_i!
export set_K_im!
export set_K_ip!
export set_K_ir!
export set_K_lp!
export set_K_lr!
export set_K_p!
export set_K_pm!
export set_K_pr!
export set_K_turb!
export set_K_ω!
export set_Ka!
export set_Kb!
export set_Kc!
export set_Kd!
export set_Kd_gov!
export set_Ke!
export set_Kf!
export set_Kg!
export set_Kh!
export set_Ki!
export set_Ki_gov!
export set_Ki_load!
export set_Ki_mw!
export set_Kl!
export set_Kp!
export set_Kp_gov!
export set_Kp_load!
export set_Ks!
export set_Kt!
export set_Kv!
export set_L_1d!
export set_L_1q!
export set_L_ad!
export set_L_aq!
export set_L_d!
export set_L_f1d!
export set_L_ff!
export set_L_q!
export set_Ld_ref!
export set_Ls_lim!
export set_PSS_flags!
export set_P_ref!
export set_R!
export set_R_1d!
export set_R_1q!
export set_R_close!
export set_R_f!
export set_R_lim!
export set_R_open!
export set_R_th!
export set_Rselect!
export set_Se!
export set_T1!
export set_T2!
export set_T3!
export set_T4!
export set_T5!
export set_T6!
export set_T7!
export set_T_AA!
export set_T_act!
export set_T_eng!
export set_T_rate!
export set_Ta!
export set_Ta_2!
export set_Ta_3!
export set_Ta_4!
export set_Ta_Tb!
export set_Tb!
export set_Tb1!
export set_Tc!
export set_Tc1!
export set_Td0_p!
export set_Td0_pp!
export set_Td_gov!
export set_Te!
export set_Tf!
export set_Tf_1!
export set_Tf_2!
export set_Tf_load!
export set_Tg!
export set_Th!
export set_Tj!
export set_Tk!
export set_Tpelec!
export set_Tq0_p!
export set_Tq0_pp!
export set_Tr!
export set_Ts!
export set_Tsa!
export set_Tsb!
export set_Tw!
export set_U0!
export set_UEL_flags!
export set_U_c!
export set_VB_max!
export set_VELM!
export set_VFE_lim!
export set_VH_max!
export set_V_lim!
export set_V_lr!
export set_V_pss!
export set_V_ref!
export set_Va_lim!
export set_Vcl!
export set_Vcu!
export set_Vdc_ref!
export set_Vf!
export set_Vi_lim!
export set_Vm_lim!
export set_Vr_lim!
export set_Wf_nl!
export set_X_th!
export set_Xd!
export set_Xd_p!
export set_Xd_pp!
export set_Xl!
export set_Xq!
export set_Xq_p!
export set_Xq_pp!
export set_Y!
export set_active_power!
export set_active_power_flow!
export set_active_power_limits!
export set_active_power_limits_from!
export set_active_power_limits_pump!
export set_active_power_limits_to!
export set_ancillary_services!
export set_angle!
export set_angle_limits!
export set_arc!
export set_area!
export set_available!
export set_b!
export set_base_power!
export set_base_voltage!
export set_battery_resistance!
export set_battery_voltage!
export set_bias!
export set_bus!
export set_bustype!
export set_cf!
export set_contributing_services!
export set_conversion_factor!
export set_db!
export set_dc_dc_inductor!
export set_dc_link_capacitance!
export set_delta_t!
export set_dynamic_injector!
export set_efficiency!
export set_energy_value!
export set_eq_p!
export set_ext!
export set_fixed!
export set_flow_limits!
export set_from!
export set_fs!
export set_fuel!
export set_fuel_flag!
export set_gate_position_limits!
export set_inflow!
export set_initial_ace!
export set_initial_energy!
export set_initial_storage!
export set_input_active_power_limits!
export set_input_code!
export set_internal_angle!
export set_internal_voltage!
export set_inv_d_fluxlink!
export set_inv_q_fluxlink!
export set_inverter_firing_angle!
export set_inverter_tap_limits!
export set_inverter_xrc!
export set_kad!
export set_kd!
export set_kffi!
export set_kffv!
export set_ki_pll!
export set_kic!
export set_kii!
export set_kiv!
export set_kp_pll!
export set_kpc!
export set_kpi!
export set_kpv!
export set_kq!
export set_kω!
export set_lf!
export set_lg!
export set_load_response!
export set_load_zone!
export set_loss!
export set_lv!
export set_magnitude!
export set_max_active_power!
export set_max_reactive_power!
export set_model!
export set_must_run!
export set_n_states!
export set_no_load!
export set_number!
export set_operation_cost!
export set_outflow!
export set_output_active_power_limits!
export set_peak_active_power!
export set_peak_reactive_power!
export set_penalty_cost!
export set_power_factor!
export set_power_trajectory!
export set_primary_shunt!
export set_prime_mover!
export set_pump_efficiency!
export set_q_nl!
export set_r!
export set_ramp_limits!
export set_ramp_limits_pump!
export set_rate!
export set_rated_current!
export set_rated_voltage!
export set_rating!
export set_rating_pump!
export set_rc_rfd!
export set_reactive_power!
export set_reactive_power_flow!
export set_reactive_power_limits!
export set_reactive_power_limits_from!
export set_reactive_power_limits_pump!
export set_reactive_power_limits_to!
export set_rectifier_firing_angle!
export set_rectifier_tap_limits!
export set_rectifier_xrc!
export set_remote_bus_control!
export set_requirement!
export set_rf!
export set_rg!
export set_rv!
export set_saturation_coeffs!
export set_services!
export set_shut_down!
export set_speed_error_signal!
export set_start_time_limits!
export set_start_types!
export set_start_up!
export set_state_of_charge_limits!
export set_states!
export set_states_types!
export set_status!
export set_storage_capacity!
export set_storage_target!
export set_switch!
export set_tap!
export set_time_at_status!
export set_time_frame!
export set_time_limits!
export set_time_limits_pump!
export set_to!
export set_valve_position_limits!
export set_variable!
export set_voltage!
export set_voltage_limits!
export set_x!
export set_α!
export set_γ_d1!
export set_γ_d2!
export set_γ_q1!
export set_γ_q2!
export set_γ_qd!
export set_γd!
export set_γq!
export set_θp!
export set_θp_rad!
export set_τ_limits!
export set_ω_lp!
export set_ωad!
export set_ωf!


text:
# code for drawing double spiral coils
"""
    double_spiral(width, height, conductor_width,
        conductor_gap, clockwise, α0, in_length, out_length)
Create a positive trace of rectangular double spiral coil of specified lengths,
conductor width, and conductor gap. The keyword arguments `clockwise` and
`α0` defines the direction of winding and the angle, respectively.
"""
function double_spiral(width, height, conductor_width,
    conductor_gap, clockwise, α0, in_length, out_length)

    minimal_length = 3 * conductor_width + 2 * conductor_gap
    if (height < minimal_length) | (width < minimal_length)
        error("""
        Cannot make a rectangular double coil with given parameters.
        Both the width and the height of the coil must be larger than
        3*conductor_width + 2*conductor_gap
        """)
    end
    # turn_direction:
    # -1 if?lum clockwise-winding coil, +1 if counter-clockwise-winding coil
    turn_direction = -2 * clockwise + 1

    point0 = -turn_direction * Point(conductor_width/2 * sin(α0),
        - conductor_width/2 * cos(α0))
    p = Path(point0, α0=α0)

    path_style = Paths.Trace(conductor_width)
    corner_style = Paths.SimpleTraceCorner()

    # first length and second length
    l1 = height - conductor_width / 2
    l2 = width - conductor_gap - 2 * conductor_width
    # create a list of lengths of straight lines
    l_list = [l1, l2]
    l = l1 - (1.5 * conductor_width + conductor_gap)
    while l >= (conductor_width + conductor_gap)
        push!(l_list, l)
        l = l_list[end-1] - 2 * (conductor_width + conductor_gap)
    end

    # length of the input line
    l_list[1] += in_length
    # scale the innermost lengths to fit nice and symmetrically
    l_list[end] = l_list[end] / 2 + (conductor_width + conductor_gap) / 2
    # middle section
    l_mid = l_list[end-1] - conductor_width - conductor_gap

    # wind in
    for n in 1:length(l_list)
        straight!(p, l_list[n], path_style)
        corner!(p, turn_direction * π/2, corner_style)
    end
    # connect (middle segment)
    if l_mid > 0 * l_mid
        straight!(p, l_mid, path_style)
    end

    l_list[1] -= in_length
    # length of the output line
    l_list[1] += out_length
    # wind out
    for n in length(l_list):-1:1
        corner!(p, -turn_direction * π/2, corner_style)
        straight!(p, l_list[n], path_style)
    end
    info("Double-spiral coil pathlength : $(pathlength(p) - in_length - out_length)")
    return p
end

"""
    double_spiral!{T,S}(cell::Cell{T,S}, width, height, conductor_width,
        conductor_gap, meta=GDSMeta(0,0); clockwise=true, α0=π/2)
Create a positive trace of rectangular double spiral coil of specified lengths,
conductor width, and conductor gap. The keyword arguments `clockwise` and
`α0` defines the direction of winding and the angle, respectively.
"""
function double_spiral!(c::Cell{T,S}, width, height, conductor_width,
    conductor_gap, meta::Meta=GDSMeta(0,0); clockwise=true, α0=π/2,
    in_length=zero(T), out_length=zero(T)) where {T<:Coordinate,S<:GDSMeta}

    p = double_spiral(width, height, conductor_width,
        conductor_gap, clockwise, α0, in_length, out_length)
    render!(c, p, meta)
    return c
end

#TODO : add method for drawing dspiral with only total length, conductor_width,
#       and conductor_gap specified


text:
#=
    This code is adapted from Flux.jl/src/functors.jl, credit to Mike Innes, et al.

    A reference to the original code can be found here:
    https://github.com/FluxML/Flux.jl/blob/v0.10.0/src/functor.jl
=#

"""
    isleaf(x)
Return true if `x` has no children according to [`parameters`](@ref).
"""
isleaf(x) = parameters(x) === ()

"""
    parameters(x)

Returns a tuple of parameters, as marked by the user during construction
"""
parameters(x) = ()

parameters(x::Tuple) = x
parameters(x::NamedTuple) = x

parameters(x::AbstractArray) = x
parameters(x::AbstractArray{<:Number}) = ()

"""
    parammap(f, x)

Maps `f` to each parameter (as specified by [`parameters`](@ref)) to the model `x`. Evaluates
    calls to [`ffunctor`](@ref).
"""
function parammap(f, x)
    func, re = ffunctor(x)
    return re(map(f, func))
end

"""
    ffunctor(x::T)

Returns a flexible functor based on the chosen [`parameters`](@ref).
"""
function ffunctor(x::T) where T
    params = parameters(x)
    re = (y) -> begin
        all_args = map(fieldnames(T)) do fn
            field = fn in params ? getfield(y, fn) : getfield(x, fn)
            return field
        end
        return constructorof(T)(all_args...)
    end
    func = (; (p => getfield(x, p) for p in params)...)
    return func, re
end
ffunctor(x::AbstractArray) = x, y -> y
ffunctor(x::AbstractArray{<:Number}) = (), _ -> x

"""
    destructure(s)

For a struct `s` which has parameters given by [`parameters`](@ref), iterates through all
    non-leaf nodes and collects the marked fields into a flat vector. This is particularly
    useful for model training or optimization, or use with `ForwardDiff`. A function which
    restructures `s` according to a flat vector is returned as the second argument.

NOTE: The flat vector representation is a `Vector{T}`. Julia will promote all entries to a
    common type if possible. This means [1, 2.0] == [1.0, 2.0] and both are Vector{Float64}.
"""
function destructure(s)
    xs = []
    fieldmap(s) do x
        x isa AbstractArray && push!(xs, x)
        x isa Number && push!(xs, [x])
        return x
    end
    return vcat(vec.(copy(xs))...), p -> restructure(s, p)
end

"""
    restructure(s, xs)

Given a struct `s` with parameters given by [`parameters`](@ref), restructures `s` according
    to the vector `xs`. In particular, `s == restructure(s, destructure(s)[1])`.
"""
function restructure(s, xs)
    i = 0
    fieldmap(s) do x
        if x isa AbstractArray
            x = reshape(xs[i.+(1:length(x))], size(x))
            i += length(x)
        elseif x isa Number
            x = xs[i+1]
            i += 1
        end
        return x
    end
end

"""
    fieldmap(f, x; exlude = isleaf)

Maps the function `f` over the fields of a `FlexibleFunctor` (equivalently, the parameter
    fields of a struct `s` given by [`parameters`](@ref)).
"""
function fieldmap(f, x; exclude = isleaf)
    y = exclude(x) ? f(x) : parammap(x -> fieldmap(f, x; exclude = exclude), x)
    return y
end


text:
# strings
x = "ppzz!#"
r = RObject(x)
@test isa(r,RObject{StrSxp})
@test length(r) == 1
@test rcopy(r) == x
@test isascii(r)
@test isascii(rcopy(r))
@test rcopy(r[1]) == x
@test isa(r[1],RObject{CharSxp})

x = "aaα₁"
r = RObject(x)
@test isa(r,RObject{StrSxp})
@test length(r) == 1
@test rcopy(r) == x
@test !isascii(rcopy(r))
@test rcopy(Symbol, r) == Symbol(x)

v = ["ap","xα⟩","pp"]
r = RObject(v)
@test isa(r,RObject{StrSxp})
@test length(r) == length(v)
@test rcopy(r) == v
@test rcopy(r[1]) == v[1]
@test rcopy(Array, r)[1] == v[1]
@test rcopy(Vector, r)[1] == v[1]
@test rcopy(Array{Compat.String}, r)[2] == v[2]
@test rcopy(Array{Symbol}, r)[2] == Symbol(v[2])
@test isa(RCall.sexp(StrSxp, :a), Ptr{StrSxp})

# logical
r = RObject(false)
@test isa(r,RObject{LglSxp})
@test length(r) == 1
@test rcopy(r) === false
@test r[1] === convert(Cint,0)

v = bitrand(10)
r = RObject(v)
@test isa(r,RObject{LglSxp})
@test length(r) == length(v)
@test size(r) == size(v)
@test rcopy(r) == v
@test rcopy(Vector, r) == v
@test rcopy(Array, r) == v
@test r[1] === convert(Cint,v[1])

# integer
x = 7
r = RObject(x)
@test isa(r,RObject{IntSxp})
@test length(r) == 1
@test size(r) == (1,)
@test rcopy(r) === convert(Cint,x)
@test r[1] === convert(Cint,x)

v = -7:3
r = RObject(v)
@test isa(r,RObject{IntSxp})
@test length(r) == length(v)
@test size(r) == size(v)
@test isa(rcopy(r), Vector{Cint})
@test rcopy(r) == collect(v)
@test rcopy(Vector, r) == collect(v)
@test r[1] === convert(Cint,v[1])
@test r[3] === convert(Cint,v[3])
r[2] = -100
@test r[2] === convert(Cint,-100)
@test r[1] === convert(Cint,v[1])
@test r[3] === convert(Cint,v[3])

m = Int[-5 2 9; 7 -8 3]
r = RObject(m)
@test isa(r,RObject{IntSxp})
@test length(r) == length(m)
@test size(r) == size(m)
@test isa(rcopy(r), Matrix{Cint})
@test rcopy(r) == m
@test rcopy(Array, r) == m
@test r[1] === convert(Cint,m[1])
@test r[3] === convert(Cint,m[3])
@test r[2,2] === convert(Cint,m[2,2])
r[2] = -100
@test r[2] === convert(Cint,-100)
r[1,3] = -101
@test r[1,3] === convert(Cint,-101)

a = rand(-20:20,2,4,5)
r = RObject(a)
@test isa(r,RObject{IntSxp})
@test length(r) == length(a)
@test size(r) == size(a)
@test isa(rcopy(r), Array{Cint,3})
@test rcopy(r) == a
@test r[1] === convert(Cint,a[1])
@test r[3] === convert(Cint,a[3])
@test r[2,3,2] === convert(Cint,a[2,3,2])

a = rand(-20:20,2,4,2,3)
r = RObject(a)
@test isa(r,RObject{IntSxp})
@test length(r) == length(a)
@test size(r) == size(a)
@test isa(rcopy(r), Array{Cint,4})
@test rcopy(r) == a
@test r[1] === convert(Cint,a[1])
@test r[3] === convert(Cint,a[3])
@test r[2,3,1,2] === convert(Cint,a[2,3,1,2])

# real
x = 7.0
r = RObject(x)
@test isa(r,RObject{RealSxp})
@test length(r) == 1
@test size(r) == (1,)
@test rcopy(r) === convert(Float64,x)
@test r[1] === convert(Float64,x)

v = -7.0:3.0
r = RObject(v)
@test isa(r,RObject{RealSxp})
@test length(r) == length(v)
@test size(r) == size(v)
@test isa(rcopy(r), Vector{Float64})
@test rcopy(r) == collect(v)
@test rcopy(Vector, r) == collect(v)
@test r[1] === convert(Float64,v[1])
@test r[3] === convert(Float64,v[3])

m = Float64[-5 2 9; 7 -8 3]
r = RObject(m)
@test isa(r,RObject{RealSxp})
@test length(r) == length(m)
@test size(r) == size(m)
@test isa(rcopy(r), Matrix{Float64})
@test rcopy(r) == m
@test rcopy(Array, r) == m
@test r[1] === convert(Float64,m[1])
@test r[3] === convert(Float64,m[3])
@test r[2,2] === convert(Float64,m[2,2])

a = rand(2,4,5)
r = RObject(a)
@test isa(r,RObject{RealSxp})
@test length(r) == length(a)
@test size(r) == size(a)
@test isa(rcopy(r), Array{Float64,3})
@test rcopy(r) == a
@test r[1] === convert(Float64,a[1])
@test r[3] === convert(Float64,a[3])
@test r[2,3,2] === convert(Float64,a[2,3,2])

a = rand(2,4,2,3)
r = RObject(a)
@test isa(r,RObject{RealSxp})
@test length(r) == length(a)
@test size(r) == size(a)
@test isa(rcopy(r), Array{Float64,4})
@test rcopy(r) == a
@test r[1] === convert(Float64,a[1])
@test r[3] === convert(Float64,a[3])
@test r[2,3,1,2] === convert(Float64,a[2,3,1,2])

# complex
x = 7.0-2.0*im
r = RObject(x)
@test isa(r,RObject{CplxSxp})
@test length(r) == 1
@test size(r) == (1,)
@test rcopy(r) === convert(Complex128,x)
@test r[1] === convert(Complex128,x)

v = randn(7)+im*randn(7)
r = RObject(v)
@test isa(r,RObject{CplxSxp})
@test length(r) == length(v)
@test size(r) == size(v)
@test isa(rcopy(r), Vector{Complex128})
@test rcopy(r) == collect(v)
@test rcopy(Vector, r) == collect(v)
@test r[1] === convert(Complex128,v[1])
@test r[3] === convert(Complex128,v[3])

m = Float64[-5 2 9; 7 -8 3] + im*Float64[-5 2 9; 7 -8 3]
r = RObject(m)
@test isa(r,RObject{CplxSxp})
@test length(r) == length(m)
@test size(r) == size(m)
@test isa(rcopy(r), Matrix{Complex128})
@test rcopy(r) == m
@test rcopy(Array, r) == m
@test r[1] === convert(Complex128,m[1])
@test r[3] === convert(Complex128,m[3])
@test r[2,2] === convert(Complex128,m[2,2])

a = rand(2,4,5)+im*randn(2,4,5)
r = RObject(a)
@test isa(r,RObject{CplxSxp})
@test length(r) == length(a)
@test size(r) == size(a)
@test isa(rcopy(r), Array{Complex128,3})
@test rcopy(r) == a
@test r[1] === convert(Complex128,a[1])
@test r[3] === convert(Complex128,a[3])
@test r[2,3,2] === convert(Complex128,a[2,3,2])

a = rand(2,4,2,3)+im*randn(2,4,2,3)
r = RObject(a)
@test isa(r,RObject{CplxSxp})
@test length(r) == length(a)
@test size(r) == size(a)
@test isa(rcopy(r), Array{Complex128,4})
@test rcopy(r) == a
@test r[1] === convert(Complex128,a[1])
@test r[3] === convert(Complex128,a[3])
@test r[2,3,1,2] === convert(Complex128,a[2,3,1,2])

# dict
d = Dict(:a=>[1, 2, 4], :b=> ["e", "d", "f"])
r = RObject(d)
@test r[:a][3] == 4
@test rcopy(r[:b][2]) == "d"
l = rcopy("list(a=1,b=c(1,3,4))")
@test l[:a] == 1
@test l[:b][3] == 4
d = RObject(Dict(1=>2))
@test Dict{Any,Any}("1" => 2) == rcopy(Dict, d)
@test Dict{Int,Int}(1=>2) == rcopy(Dict{Int,Int}, d)

# function
function funk(x,y)
    x+y
end
f1 = RObject(funk)
@test rcopy(Function, f1)(1,2) == 3
@test rcopy(Function, f1.p)(1,2) == 3

# misc
a = RObject(rand(10))
@test length(rcopy(Any, a)) == 10
@test typeof(RCall.sexp(Cint, 1)) == Cint
@test typeof(RCall.sexp(Float64, 1)) == Float64
@test typeof(RCall.sexp(Complex128, 1)) == Complex128
@test typeof(rcopy(Vector{Float64}, a.p)) == Vector{Float64}
b = RObject(true)
@test rcopy(Int32(1)) == 1
@test rcopy(Cint, Int32(1)) == 1
@test rcopy(Cint, b.p) == 1
@test rcopy(Vector{Cint}, b.p) == [1]
@test rcopy(Array{Cint}, b.p) == [1]
@test rcopy(Vector{Bool}, b.p) == [true]
@test rcopy(BitVector, b.p) == [true]

#RCall.rlang_formula(parse("a+b"))
@test RCall.rlang_formula(:a) == :a


text:
using Test

@testset "StatsFunctions.jl" begin
    include("StatsFunctionsTest.jl")
end

@testset "Backward.jl" begin
    include("BackwardTest.jl")
end

@testset "Hierarchical Gaussian Gradient" begin
    include("HierGaussLikelihood.jl")
end


text:
function rsi{T}(ta::TimeArray{T,1}, n::Int; wilder=false)

    ret = [0; diff(ta.values)]
    ups = zeros(length(ta))
    dns = zeros(length(ta))

    for i=1:length(ta)
        if ret[i] >= 0
            ups[i] += ret[i]
        else
            dns[i] += ret[i]
        end
     end

    if  wilder
        upsema = ema(ups, n, wilder=true)
        dnsema = abs(ema(dns, n, wilder=true))
        rs     = upsema ./ dnsema

    else
        upsema  = ema(ups, n)
        dnsema  = abs(ema(dns, n))
        rs      = upsema ./ dnsema
    end

    res  = 100 .- (100./(1 .+ rs))

    TimeArray(ta.timestamp[n:end], res, ["rsi"], ta.meta)
end

rsi{T}(ta::TimeArray{T,1}) = rsi(ta, 14)

function macd{T}(ta::TimeArray{T,1}, fast::Int, slow::Int, signal::Int)
    mcd = ema(ta, fast) .- ema(ta, slow)
    sig = ema(mcd, signal)
    merge(mcd, sig, col_names=["macd", "signal"])
end

macd{T}(ta::TimeArray{T,1}) = macd(ta, 12, 26, 9)

# function cci{T,N,M}(ohlc::TimeArray{T,N,M}, ma::Int, c::Float64)
#   	typ     = typical(ohlc)
#     sma_typ = sma(typ, ma)
#     diff    = typ .- sma_typ
#     #mead   = mad1(typ)
#     mead    = moving(typ, mad1, 20)
#     divisor = mead .* c
#     vals    = diff ./ divisor
#   	TimeArray(diff.timestamp, vals.values, ["cci"])
# end
#
# cci{T,N,M}(ohlc::TimeArray{T,N,M}) = cci(ohlc, 20, .015)


text:
@testset "Test Fractional Integral Equations SpectralUltraspherical method" begin
    tspan = LinRange(-1, 1, 5)
    prob = FIEProblem([1, 1], [1, 0.5], 1, tspan)
    sol = solve(prob, 20, SpectralUltraspherical())

    @test isapprox(sol.u, [ 1.0
    0.5231565837302468
    0.42758357615580733
    0.37316567427801584
    0.3362040024463395]; atol=1e-3)
end


text:
"""
given `iso_vals` and iso, return an 8 bit value corresping
to each corner of a cube. 0 indicates in the isosurface
and 1 indicates outside the surface
"""
@inline function _get_cubeindex(iso_vals, iso)
    cubeindex = iso_vals[1] < iso ? 0x01 : 0x00
    iso_vals[2] < iso && (cubeindex |= 0x02)
    iso_vals[3] < iso && (cubeindex |= 0x04)
    iso_vals[4] < iso && (cubeindex |= 0x08)
    iso_vals[5] < iso && (cubeindex |= 0x10)
    iso_vals[6] < iso && (cubeindex |= 0x20)
    iso_vals[7] < iso && (cubeindex |= 0x40)
    iso_vals[8] < iso && (cubeindex |= 0x80)
    cubeindex
end


text:
@benchset "abstract_gp" begin
    for N in [10,]
        f = GP(CustomMean(sin), SEKernel(), GPC())
        x = randn(N)
        fx = f(x)
        y = rand(fx)
        @benchset "N=$N" begin
            create_benchmarks("logpdf SEKernel $N (x, y)", true, (x, y)->logpdf(f(x), y), x, y)
            create_benchmarks("rand SEKernel $N", true, x->rand(f(x)), x)
            create_benchmarks("cov(f(x)) $N", true, x->cov(f(x)), x)
        end
    end
end


text:
__precompile__(true)

module BayesianNonparametrics

  using Distributions#, Combinatorics,
  using Clustering, ProgressMeter, StatsFuns

  # include Base modules
  using SparseArrays, Distributed, LinearAlgebra, Random

  include("common.jl")
  include("math.jl")
  include("utils.jl")
  include("datasets.jl")
  include("distributions.jl")
  include("distfunctions.jl")
  include("inits.jl")
  include("inference.jl")
  include("models.jl")
  include("dpmm.jl")
  include("hdp.jl")
  include("vipointestimate.jl")

end # module


text:
# by lazarusA # HIDE
#using GLMakie, Random # HIDE
using CairoMakie, Random
CairoMakie.activate!() #HIDE
let
    Random.seed!(123)
    fig = Figure(resolution = (400, 427))
    ax = Axis(fig,  aspect = 1, xlabel = "x", ylabel = "y")
    hmap = heatmap!(rand(20,20), colormap = :CMRmap)
    cbar = Colorbar(fig, hmap, label = "values", height = 15, vertical = false,
                ticksize=15, tickalign = 1, width = Relative(3.55/4))
    fig[1, 1] = cbar
    fig[2, 1] = ax
    rowgap!(fig.layout, 7)
    fig
    save(joinpath(@__DIR__, "output", "heatmapCbarTop.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["CairoMakie", "Random"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

using DataFrames
using Datasets
using RLESUtils, PGFPlotUtils
using TikzPictures

const COLNAMEMAP = [:id=>:id, :A=>:A, :L=>:L, :_chi=>:chi, :_beta=>:beta,
    :C_1=>:C1, :C_2=>:C2,:v_1=>:v1, :v_2=>:v2, :_dot_v_1=>:v1d, :_dot_v_2=>:v2d,
    :_dot_h_1=>:h1d, :_dot_h_2=>:h2d, :_dot_psi_1=>:psi1d, :_dot_psi_2=>:psi2d,
    :hmd=>:hmd, :vmd=>:vmd]

function initial_to_DataFrame(init_file::AbstractString)
    D = readtable(init_file; separator=' ')
    for (from,to) in COLNAMEMAP
        if from != to
            rename!(D, from, to)
        end
    end
    D
end

"""
Repairs the id number so that it matches the row number.  This is useful for fixing
the id numbers after, for example, manually replacing the first encounter with the
standard header.
"""
function repair_ids(init_file::AbstractString; skipfirst::Bool=true, delim=' ')
    toks = splitext(init_file)
    outfile = toks[1] * "_repaired" * toks[2]
    fin = open(init_file, "r")
    fout = open(outfile, "w")
    i = 1
    if skipfirst #header
        line = readline(fin)
        print(fout, line)
    end
    for line in readlines(fin)
        toks = split(line, delim)
        toks[1] = "$i"
        join(fout, toks, delim)
        i += 1
    end
    close(fin)
    close(fout)
end

function llcem_initial_study(init_file::AbstractString, data_name::AbstractString;
    outfileroot::AbstractString="llcem_initial", format::Symbol=:TEXPDF)
    m = load_meta(data_name)
    nmac_ids = find(m[:nmac])
    init = initial_to_DataFrame(init_file)
    td = histogram_all_cols_sbs(init[nmac_ids,:], init;
        discretization=:sqrt,
        datanames=["nmacs", "all"])
    plot_tikz(outfileroot, td, format)
end

max_stars_repo_name: sisl/RLESCAS

max_issues_repo_name: sisl/RLESCAS

max_forks_repo_name: sisl/RLESCAS


text:
# From https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/naif_ids.html
# TODO do we need the rest?
export naif_body_codes, naif_reference_codes
naif_body_codes = Dict(
                       # Barycenters
                       "SOLAR_SYSTEM_BARYCENTER" => 0, "SSB" => 0,
                       "SOLAR SYSTEM BARYCENTER" => 0, "MERCURY_BARYCENTER" => 1,
                       "MERCURY BARYCENTER" => 1, "VENUS_BARYCENTER" => 2,
                       "VENUS BARYCENTER" => 2, "EARTH_BARYCENTER" => 3, "EMB" => 3,
                       "EARTH MOON BARYCENTER" => 3, "EARTH-MOON BARYCENTER" => 3,
                       "EARTH BARYCENTER" => 3, "MARS_BARYCENTER" => 4,
                       "MARS BARYCENTER" => 4, "JUPITER_BARYCENTER" => 5,
                       "JUPITER BARYCENTER" => 5, "SATURN_BARYCENTER" => 6,
                       "SATURN BARYCENTER" => 6, "URANUS_BARYCENTER" => 7,
                       "URANUS BARYCENTER" => 7, "NEPTUNE_BARYCENTER" => 8,
                       "NEPTUNE BARYCENTER" => 8, "PLUTO_BARYCENTER" => 9,
                       "PLUTO BARYCENTER" => 9, "SUN" => 10,
                       # Planets and Satellites
                       "MERCURY" => 199, "VENUS" => 299, "EARTH" => 399, "MOON" => 301,
                       "MARS" => 499, "PHOBOS" => 401, "DEIMOS" => 402, "JUPITER" => 599,
                       "IO" => 501, "EUROPA" => 502, "GANYMEDE" => 503, "CALLISTO" => 504,
                       "AMALTHEA" => 505, "HIMALIA" => 506, "ELARA" => 507,
                       "PASIPHAE" => 508, "SINOPE" => 509, "LYSITHEA" => 510,
                       "CARME" => 511, "ANANKE" => 512, "LEDA" => 513, "THEBE" => 514,
                       "ADRASTEA" => 515, "METIS" => 516, "CALLIRRHOE" => 517,
                       "THEMISTO" => 518, "MAGACLITE" => 519, "TAYGETE" => 520,
                       "CHALDENE" => 521, "HARPALYKE" => 522, "KALYKE" => 523,
                       "IOCASTE" => 524, "ERINOME" => 525, "ISONOE" => 526,
                       "PRAXIDIKE" => 527, "AUTONOE" => 528, "THYONE" => 529,
                       "HERMIPPE" => 530, "AITNE" => 531, "EURYDOME" => 532,
                       "EUANTHE" => 533, "EUPORIE" => 534, "ORTHOSIE" => 535,
                       "SPONDE" => 536, "KALE" => 537, "PASITHEE" => 538, "HEGEMONE" => 539,
                       "MNEME" => 540, "AOEDE" => 541, "THELXINOE" => 542, "ARCHE" => 543,
                       "KALLICHORE" => 544, "HELIKE" => 545, "CARPO" => 546,
                       "EUKELADE" => 547, "CYLLENE" => 548, "KORE" => 549, "HERSE" => 550,
                       "DIA" => 553, "SATURN" => 699, "MIMAS" => 601, "ENCELADUS" => 602,
                       "TETHYS" => 603, "DIONE" => 604, "RHEA" => 605, "TITAN" => 606,
                       "HYPERION" => 607, "IAPETUS" => 608, "PHOEBE" => 609, "JANUS" => 610,
                       "EPIMETHEUS" => 611, "HELENE" => 612, "TELESTO" => 613,
                       "CALYPSO" => 614, "ATLAS" => 615, "PROMETHEUS" => 616,
                       "PANDORA" => 617, "PAN" => 618, "YMIR" => 619, "PAALIAQ" => 620,
                       "TARVOS" => 621, "IJIRAQ" => 622, "SUTTUNGR" => 623, "KIVIUQ" => 624,
                       "MUNDILFARI" => 625, "ALBIORIX" => 626, "SKATHI" => 627,
                       "ERRIAPUS" => 628, "SIARNAQ" => 629, "THRYMR" => 630, "NARVI" => 631,
                       "METHONE" => 632, "PALLENE" => 633, "POLYDEUCES" => 634,
                       "DAPHNIS" => 635, "AEGIR" => 636, "BEBHIONN" => 637,
                       "BERGELMIR" => 638, "BESTLA" => 639, "FARBAUTI" => 640,
                       "FENRIR" => 641, "FORNJOT" => 642, "HATI" => 643, "HYRROKKIN" => 644,
                       "KARI" => 645, "LOGE" => 646, "SKOLL" => 647, "SURTUR" => 648,
                       "ANTHE" => 649, "JARNSAXA" => 650, "GREIP" => 651, "TARQEQ" => 652,
                       "AEGAEON" => 653, "URANUS" => 799, "ARIEL" => 701, "UMBRIEL" => 702,
                       "TITANIA" => 703, "OBERON" => 704, "MIRANDA" => 705,
                       "CORDELIA" => 706, "OPHELIA" => 707, "BIANCA" => 708,
                       "CRESSIDA" => 709, "DESDEMONA" => 710, "JULIET" => 711,
                       "PORTIA" => 712, "ROSALIND" => 713, "BELINDA" => 714, "PUCK" => 715,
                       "CALIBAN" => 716, "SYCORAX" => 717, "PROSPERO" => 718,
                       "SETEBOS" => 719, "STEPHANO" => 720, "TRINCULO" => 721,
                       "FRANCISCO" => 722, "MARGARET" => 723, "FERDINAND" => 724,
                       "PERDITA" => 725, "MAB" => 726, "CUPID" => 727, "NEPTUNE" => 899,
                       "TRITON" => 801, "NEREID" => 802, "NAIAD" => 803, "THALASSA" => 804,
                       "DESPINA" => 805, "GALATEA" => 806, "LARISSA" => 807,
                       "PROTEUS" => 808, "HALIMEDE" => 809, "PSAMATHE" => 810, "SAO" => 811,
                       "LAOMEDEIA" => 812, "NESO" => 813, "PLUTO" => 999, "CHARON" => 901,
                       "NIX" => 902, "HYDRA" => 903, "KERBEROS" => 904, "STYX" => 905)

# https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/frames.html
naif_reference_codes = Dict("J2000" => 1, "B1950" => 2, "FK4" => 3, "DE-118" => 4,
                            "DE-96" => 5, "DE-102" => 6, "DE-108" => 7, "DE-111" => 8,
                            "DE-114" => 9, "DE-122" => 10, "DE-125" => 11, "DE-130" => 12,
                            "GALACTIC" => 13, "DE-200" => 14, "DE-202" => 15,
                            "MARSIAU" => 16, "ECLIPJ2000" => 17, "ECLIPB1950" => 18,
                            "DE-140" => 19, "DE-142" => 20, "DE-143" => 21)


text:
# ---
# title: Custom median filters
# cover: assets/median.gif
# author: Johnny Chen
# date: 2020-09-23
# ---

# With median filter as an example, this demo shows how you could construct a custom kernel and
# pass it to [`mapwindow`](@ref) for common stencil operations.

# Unlike other image processing toolboxes, `ImageFiltering.jl` does not provide you an out-of-box
# function for median filters. It instead provides the more general [`mapwindow`](@ref) function,
# which allows you to apply any function to the window around each pixel.

using ImageFiltering, ImageCore, ImageShow # or you could just `using Images`
using TestImages
using Statistics
using Random #hide
Random.seed!(0) #hide

img = Float64[isodd(i) - isodd(j) for i = 1:5, j = 1:5]
img[3, 3] = 1000
img #hide #md

#-

patch_size = (3, 3)
imgm = mapwindow(median, img, patch_size)

# `mapwindow` provides a high-level interface to loop over the image and call a function (`median`
# in this demo) on each patch/window. It's performing an operation that is nearly equivalent to the
# loop below, albeit more efficiently:

## For simplicity, border condition is not included here.
imgm = zeros(axes(img))
R = CartesianIndices(img)
I_first, I_last = first(R), last(R)
Δ = CartesianIndex(patch_size .÷ 2)
for I in R
    patch = max(I_first, I-Δ):min(I_last, I+Δ)
    imgm[I] = median(img[patch])
end
imgm

# Compared to this hand-written loop, `mapwindow` offers additional flexibility in handling the
# borders and allows you to process only a subset of the windows. For some functions (e.g., `min`,
# `max`, and `extrema`), `mapwindow` uses a custom implementation that is far more efficient than
# naively computing the function "freshly" on each window.

# When the input array has `NaN` or some unwanted pixel values, a pre-filtering process is needed to
# exclude them first. This can be done quite easily by compositing a given kernel operation and a
# `filter` operation. Let's still take `median` filter as an example.

img[2,2] = NaN
imgm = mapwindow(median, img, (3,3))

# `NaN` has polluted the output, which is usually not wanted. This can be fixed quite easily by
# compositing a new filter kernel. This way, we only compute the median value on the non-NaN subset
# of each patch/window.

_median(patch) = median(filter(x->!isnan(x), patch))
imgm = mapwindow(_median, img, patch_size)

# ## Impulse (Salt and Pepper) noise removal

# Median filters are quite robust to outliers (unusual values), this property can be used to remove
# salt&pepper noise. An image with `n`-level salt&papper noise is defined as: each pixel `p` has
# `n/2` probabilty to be filled by `0`(the minimal gamut value), `n/2` probabilty to be filled by
# `1`(the maximal gamut value), and `1-n` probablity to be undamaged.

img = testimage("cameraman")

n = 0.2 # noise level
noisy_img = map(img) do p
    sp = rand()
    if sp < n/2
        eltype(img)(gamutmin(eltype(img))...)
    elseif sp < n
        eltype(img)(gamutmax(eltype(img))...)
    else
        p
    end
end

denoised_img = mapwindow(median!, noisy_img, (5, 5))

mosaicview(img, noisy_img, denoised_img; nrow=1)

## save covers #src
using ImageMagick #src
mkpath("assets") #src
ImageMagick.save("assets/median.gif", cat(noisy_img, denoised_img; dims=3); fps=1) #src


text:
module FlightSims

using Reexport
@reexport using FSimBase
using OrdinaryDiffEq: Tsit5  # default solver
@reexport using FSimZoo

using UnPack

export ned2enu, enu2ned

include("APIs/APIs.jl")
include("utils/utils.jl")

end


text:
# variational autoencoder with Gaussian observed and latent variables
mutable struct VAE{T}
    # encoder / recognizer
    We1::AbstractMatrix{T}  # encoder: layer 1 weight
    be1::AbstractVector{T}  # encoder: layer 1 bias
    We2::AbstractMatrix{T}  # encoder: layer 2 weight
    be2::AbstractVector{T}  # encoder: layer 2 bias
    We3::AbstractMatrix{T}  # encoder: layer 3 mu weight
    be3::AbstractVector{T}  # encoder: layer 3 mu bias
    We4::AbstractMatrix{T}  # encoder: layer 3 log(sigma^2) weight
    be4::AbstractVector{T}  # encoder: layer 3 log(sigma^2) bias
    # decoder / generator
    Wd1::AbstractMatrix{T}  # decoder: layer 1 weight
    bd1::AbstractVector{T}  # decoder: layer 1 bias
    Wd2::AbstractMatrix{T}  # decoder: layer 2 weight
    bd2::AbstractVector{T}  # decoder: layer 2 bias
    Wd3::AbstractMatrix{T}  # decoder: layer 3 weight
    bd3::AbstractVector{T}  # decoder: layer 3 bias
end

function Base.show(io::IO, m::VAE{T}) where T
    print(io, "VAE{$T}($(size(m.We1,2)), $(size(m.We1,1)), $(size(m.We2,1)), " *
          "$(size(m.We3,1)), $(size(m.Wd1,1)), $(size(m.Wd2,1)), $(size(m.Wd3,1)))")
end

VAE{T}(n_inp, n_he1, n_he2, n_z, n_hd1, n_hd2, n_out) where T =
    VAE{T}(
        # encoder
        xavier_init(n_he1, n_inp),
        zeros(n_he1),
        xavier_init(n_he2, n_he1),
        zeros(n_he2),
        xavier_init(n_z, n_he2),
        zeros(n_z),
        xavier_init(n_z, n_he2),
        zeros(n_z),
        # decoder
        xavier_init(n_hd1, n_z),
        zeros(n_hd1),
        xavier_init(n_hd2, n_hd1),
        zeros(n_hd2),
        xavier_init(n_out, n_hd2),
        zeros(n_out)
    )

# empty constructor is required by XGrad
VAE{T}() where T = VAE{T}(1, 1, 1, 1, 1, 1, 1)

function encode(m::VAE, x)
    he1 = softplus.(m.We1 * x .+ m.be1)
    he2 = softplus.(m.We2 * he1 .+ m.be2)
    mu = m.We3 * he2 .+ m.be3
    log_sigma2 = m.We4 * he2 .+ m.be4
    return mu, log_sigma2
end

function decode(m::VAE, z)
    hd1 = softplus.(m.Wd1 * z .+ m.bd1)
    hd2 = softplus.(m.Wd2 * hd1 .+ m.bd2)
    x_rec = logistic.(m.Wd3 * hd2 .+ m.bd3)
    return x_rec
end

# TODO: make cost with Gaussian decoder
# TODO: figure out why KLD is different from the one in papers
function vae_cost(m::VAE, eps, x)
    mu, log_sigma2 = encode(m, x)
    z = mu .+ sqrt.(exp.(log_sigma2)) .* eps
    x_rec = decode(m, z)
    # loss
    rec_loss = -sum(x .* log.(1e-10 + x_rec) + (1 - x) .* log.(1e-10 + 1 - x_rec), 1)
    KLD = -0.5 * sum(1 + log_sigma2 .- mu .^ 2 - exp.(log_sigma2), 1)
    cost = mean(rec_loss .+ KLD)
end