text:
"""
	characteristicMatrix( FV::Cells )::ChainOp

Binary matrix representing by rows the `p`-cells of a cellular complex.
The input parameter must be of `Cells` type. Return a sparse binary matrix,
providing the basis of a ``Chain`` space of given dimension. Notice that the
number of columns is equal to the number of vertices (0-cells).

# Example

```julia
V,(VV,EV,FV,CV) = cuboid([1.,1.,1.], true);

julia> Matrix(characteristicMatrix(FV))
6×8 Array{Int8,2}:
1  1  1  1  0  0  0  0
0  0  0  0  1  1  1  1
1  1  0  0  1  1  0  0
0  0  1  1  0  0  1  1
1  0  1  0  1  0  1  0
0  1  0  1  0  1  0  1

julia> Matrix(characteristicMatrix(CV))
1×8 Array{Int8,2}:
1  1  1  1  1  1  1  1

julia> Matrix(characteristicMatrix(EV))
12×8 Array{Int8,2}:
1  1  0  0  0  0  0  0
0  0  1  1  0  0  0  0
0  0  0  0  1  1  0  0
0  0  0  0  0  0  1  1
1  0  1  0  0  0  0  0
0  1  0  1  0  0  0  0
0  0  0  0  1  0  1  0
0  0  0  0  0  1  0  1
1  0  0  0  1  0  0  0
0  1  0  0  0  1  0  0
0  0  1  0  0  0  1  0
0  0  0  1  0  0  0  1
```
"""
function characteristicMatrix( FV::Cells )::ChainOp
	I,J,V = Int64[],Int64[],Int8[]
	for f=1:length(FV)
		for k in FV[f]
		push!(I,f)
		push!(J,k)
		push!(V,1)
		end
	end
	M_2 = sparse(I,J,V)
	return M_2
end

"""
	boundary_1( EV::Cells )::ChainOp

Computation of sparse signed boundary operator ``C_1 -> C_0``.

# Example
```julia
julia> V,(VV,EV,FV,CV) = cuboid([1.,1.,1.], true);

julia> EV
12-element Array{Array{Int64,1},1}:
[1, 2]
[3, 4]
...
[2, 6]
[3, 7]
[4, 8]

julia> boundary_1( EV::Cells )
8×12 SparseMatrixCSC{Int8,Int64} with 24 stored entries:
[1 ,  1]  =  -1
[2 ,  1]  =  1
[3 ,  2]  =  -1
...       ...
[7 , 11]  =  1
[4 , 12]  =  -1
[8 , 12]  =  1

julia> Matrix(boundary_1(EV::Cells))
8×12 Array{Int8,2}:
-1   0   0   0  -1   0   0   0  -1   0   0   0
1   0   0   0   0  -1   0   0   0  -1   0   0
0  -1   0   0   1   0   0   0   0   0  -1   0
0   1   0   0   0   1   0   0   0   0   0  -1
0   0  -1   0   0   0  -1   0   1   0   0   0
0   0   1   0   0   0   0  -1   0   1   0   0
0   0   0  -1   0   0   1   0   0   0   1   0
0   0   0   1   0   0   0   1   0   0   0   1
```
"""
function boundary_1( EV::Cells )::ChainOp
	out = characteristicMatrix(EV)'
	for e = 1:length(EV)
		out[EV[e][1],e] = -1
	end
	return out
end

"""
	coboundary_0(EV::Cells)

Return the `coboundary_0` signed operator `C_0` -> `C_1`.
"""
coboundary_0(EV::Cells) = convert(ChainOp,transpose(boundary_1(EV::Cells)))

"""
	fix_redundancy(target_mat, cscFV,cscEV)

*Fix the coboundary_1 matrix*, generated by sparse matrix product ``FV * EV^t``, for complexes with some *non-convex cells*. This approach can be used when both `EV` and `FV` of the cellular complex are known. It is exact when cells are convex. Maybe non-exact, introducing spurious incidence coefficients (``redundancies``), when adjacent faces share an edge combinatorially, but not geometrically. This happen when an edge is on the boundary of face A, but only its vertices are on the boundary of face B.  TODO: Similar situations may appear when computing algebraically CF as product of known CV and FV, with non-convex cells.

In order to remove such ``redundancies``, the Euler characteristic of 2-sphere is used, where V-E+F=2. Since we have F=2 (inner and outer face), ``V=E`` must hold, and `d=E-V` is the (non-negative) ``defect`` number, called `nfixs` in the code. It equates the number of columns `edges`
whose sum is greater than 2 for the considered row (face). Remember the in a ``d``-complex, *including* the ``outer cell``, all ``(d-1)``-faces must be shared by exactly 2 ``d``-faces. Note that `FV` *must include* the row of outer shell (exterior face).

# Example

```julia
FV = [[1,2,3,4,5,17,16,12],
[1,2,3,4,6,7,8,9,10,11,12,13,14,15],
[4,5,9,11,12,13,14,15,16,17],
[2,3,6,7], [8,9,10,11]]

FE = [[1,2,3,4,9,20,17,5],
[1,6,10,7,3,8,11,12,14,15,19,18,16,5],
[4,9,20,17,16,18,19,15,13,8],
[2,10,6,7], [11,12,13,14]]

EV = [[1,2],[2,3],[3,4],[4,5],[1,12],[2,6],[3,7],[4,9],[5,17],[6,7],[8,9],
[8,10],[9,11],[10,11],[11,15],[12,13],[12,16],[13,14],[14,15],[16,17]]

V = [0   2   5   7  10   2   5   3   7  3  7  0  3  3  7  0  10;
    16  16  16  16  16  13  13  11  11  8  8  5  5  2  2  0   0]

cscFE = u_coboundary_1( FV::Cells, EV::Cells, false);
Matrix(cscFE)
```

Notice that there are two columns (2 and 13) with 3 ones, hence (3-2)+(3-2)=2 defects to fix. The fixed complex can be shown graphically as:

```julia

VV = [[k] for k in 1:size(V,2)];
using Plasm
Plasm.view( Plasm.numbering(3)((V,[VV, EV, FV])) )
```
"""
function fix_redundancy(target_mat, cscFV,cscEV) # incidence numbers > 2#E
	nfixs = 0
	faces2fix = []
	edges2fix = []
	# target_mat and cscFV (ref_mat) should have same length per row !
	for face = 1:size(target_mat,1)
		nedges = sum(findnz(target_mat[face,:])[2])
		nverts = sum(findnz(cscFV[face,:])[2])
		if nedges != nverts
			nfixs += nedges - nverts
			#println("face $face, nedges=$nedges, nverts=$nverts")
			push!(faces2fix,face)
		end
	end
	for edge = 1:size(target_mat,2)
		nfaces = sum(findnz(target_mat[:,edge])[2])
		if nfaces > 2
			#println("edge $edge, nfaces=$nfaces")
			push!(edges2fix,edge)
		end
	end
	#println("nfixs=$nfixs")
	pairs2fix = []
	for fh in faces2fix		# for each face to fix
		for ek in edges2fix		# for each edge to fix
			if target_mat[fh, ek]==1	# edge to fix \in face to fix
				v1,v2 = findnz(cscEV[ek,:])[1]
				weight(v) = length( intersect(
							findnz(cscEV[:,v])[1], findnz(target_mat[fh,:])[1] ))
				if weight(v1)>2 && weight(v2)>2
					#println("(fh,ek) = $((fh,ek))")
					push!( pairs2fix, (fh,ek) )
				end
			end
		end
	end
	for (fh,ek) in pairs2fix
		target_mat[fh, ek] = 0
	end
	cscFE = dropzeros(target_mat)
	@assert nnz(cscFE) == 2*size(cscFE,2)
	return cscFE
end
function fix_lack(target_mat, cscFV,cscEV) # incidence numbers < 2#E
end

"""
	u_coboundary_1( FV::Cells, EV::Cells, convex=true)::ChainOp

Compute the sparse *unsigned* coboundary_1 operator ``C_1 -> C_2``.
Notice that the output matrix is `m x n`, where `m` is the number of faces, and `n`
is the number of edges.

# Examples

##  Cellular complex with convex-cells, and without outer cell

```julia
julia> V,(VV,EV,FV,CV) = cuboid([1.,1.,1.], true);

julia> u_coboundary_1(FV,EV)
6×12 SparseMatrixCSC{Int8,Int64} with 24 stored entries:
[1 ,  1]  =  1
[3 ,  1]  =  1
[1 ,  2]  =  1
[4 ,  2]  =  1
...		...
[4 , 11]  =  1
[5 , 11]  =  1
[4 , 12]  =  1
[6 , 12]  =  1

julia> Matrix(u_coboundary_1(FV,EV))
6×12 Array{Int8,2}:
1  1  0  0  1  1  0  0  0  0  0  0
0  0  1  1  0  0  1  1  0  0  0  0
1  0  1  0  0  0  0  0  1  1  0  0
0  1  0  1  0  0  0  0  0  0  1  1
0  0  0  0  1  0  1  0  1  0  1  0
0  0  0  0  0  1  0  1  0  1  0  1

julia> unsigned_boundary_2 = u_coboundary_1(FV,EV)';
```

Compute the *Unsigned* `coboundary_1` operator matrix as product of two
sparse characteristic matrices.

##  Cellular complex with non-convex cells, and with outer cell

```julia
FV = [[1,2,3,4,5,17,16,12], # outer cell
[1,2,3,4,6,7,8,9,10,11,12,13,14,15],
[4,5,9,11,12,13,14,15,16,17],
[2,3,6,7], [8,9,10,11]]

EV = [[1,2],[2,3],[3,4],[4,5],[1,12],[2,6],[3,7],[4,9],[5,17],[6,7],[8,9],
[8,10],[9,11],[10,11],[11,15],[12,13],[12,16],[13,14],[14,15],[16,17]]

out = u_coboundary_1( FV::Cells, EV::Cells, false)
```
In case of expected 2-chains with non-convex cells, instance the method with
`convex = false`, in order to fix a possible redundancy of incidence values, induced by computation through multiplication of characteristic matrices. (Look at columns
2 and 13 before, generated by default).
"""
function u_coboundary_1( FV::Cells, EV::Cells, convex=true::Bool)::ChainOp
	cscFV = characteristicMatrix(FV)
	cscEV = characteristicMatrix(EV)
	out = u_coboundary_1( cscFV::ChainOp, cscEV::ChainOp, convex::Bool)
	return out
end
function u_coboundary_1( cscFV::ChainOp, cscEV::ChainOp, convex=true::Bool)::ChainOp
	temp = cscFV * cscEV'
	I,J,Val = Int64[],Int64[],Int8[]
	for j=1:size(temp,2)
		for i=1:size(temp,1)
			if temp[i,j] == 2
				push!(I,i)
				push!(J,j)
				push!(Val,1)
			end
		end
	end
	cscFE = SparseArrays.sparse(I,J,Val)
	if !convex
		cscFE = fix_redundancy(cscFE,cscFV,cscEV)
	end
	return cscFE
end

"""
	coboundary_1( FV::Cells, EV::Cells)::ChainOp

Generate the *unsigned* sparse matrix of the coboundary_1 operator.
For each row, start with the first incidence number positive (i.e. assign the orientation of the first edge to the 1-cycle of the face), then bounce back and forth between vertex columns/rows of EV and FE.

# Example

copFE = coboundary_1(FV::Cells, EV::Cells)

julia> Matrix(cscFE)
5×20 Array{Int8,2}:
 1  1  1  1  1  0  0  0  1  0  0  0  0  0  0  0  1  0  0  1
 1  0  1  0  1  1  1  1  0  1  1  1  0  1  1  1  0  1  1  0
 0  0  0  1  0  0  0  1  1  0  0  0  1  0  1  1  1  1  1  1
 0  1  0  0  0  1  1  0  0  1  0  0  0  0  0  0  0  0  0  0
 0  0  0  0  0  0  0  0  0  0  1  1  1  1  0  0  0  0  0  0
"""
function coboundary_1(FV::Array{Array{Int64,1},1}, EV::Array{Array{Int64,1},1}) # (::Cells, ::Cells)
	copFV = lar2cop(FV)
	I,J,Val = findnz(Lar.lar2cop(EV))
	copVE = sparse(J,I,Val)
	triples = hcat([[i,j,1]  for (i,j,v)  in zip(findnz(copFV * copVE)...) if v==2]...)
	I,J,Val = triples[1,:], triples[2,:], triples[3,:]
	Val = convert(Array{Int8,1},Val)
	copFE = sparse(I,J,Val)
	return copFE
end

function coboundary_1( V::Points, FV::Cells, EV::Cells, convex=true::Bool, exterior=false::Bool)::ChainOp
	# generate unsigned operator's sparse matrix
	cscFV = characteristicMatrix(FV)
	cscEV = characteristicMatrix(EV)
	##if size(V,1) == 3
		##copFE = u_coboundary_1( FV::Cells, EV::Cells )
	##elseif size(V,1) == 2
		# greedy generation of incidence number signs
		copFE = coboundary_1( V, cscFV::ChainOp, cscEV::ChainOp, convex, exterior)
	##end
	return copFE
end

function coboundary_1( V::Points, cscFV::ChainOp, cscEV::ChainOp, convex=true::Bool, exterior=false::Bool )::ChainOp

	cscFE = u_coboundary_1( cscFV::ChainOp, cscEV::ChainOp, convex)
	EV = [findnz(cscEV[k,:])[1] for k=1:size(cscEV,1)]
	cscEV = sparse(coboundary_0( EV::Cells ))
	for f=1:size(cscFE,1)
		chain = findnz(cscFE[f,:])[1]	#	dense
		cycle = spzeros(Int8,cscFE.n)	#	sparse

		edge = findnz(cscFE[f,:])[1][1]; sign = 1
		cycle[edge] = sign
		chain = setdiff( chain, edge )
		while chain != []
			boundary = sparse(cycle') * cscEV
			_,vs,vals = findnz(dropzeros(boundary))

			rindex = vals[1]==1 ? vf = vs[1] : vf = vs[2]
			r_boundary = spzeros(Int8,cscEV.n)	#	sparse
			r_boundary[rindex] = 1
			r_coboundary = cscEV * r_boundary
			r_edge = intersect(findnz(r_coboundary)[1],chain)[1]
			r_coboundary = spzeros(Int8,cscEV.m)	#	sparse
			r_coboundary[r_edge] = EV[r_edge][1]<EV[r_edge][2] ? 1 : -1

			lindex = vals[1]==-1 ? vi = vs[1] : vi = vs[2]
			l_boundary = spzeros(Int8,cscEV.n)	#	sparse
			l_boundary[lindex] = -1
			l_coboundary = cscEV * l_boundary
			l_edge = intersect(findnz(l_coboundary)[1],chain)[1]
			l_coboundary = spzeros(Int8,cscEV.m)	#	sparse
			l_coboundary[l_edge] = EV[l_edge][1]<EV[l_edge][2] ? -1 : 1

			if r_coboundary != -l_coboundary  # false iff last edge
				# add edge to cycle from both sides
				rsign = rindex == EV[r_edge][1] ? 1 : -1
				lsign = lindex == EV[l_edge][2] ? -1 : 1
				cycle = cycle + rsign * r_coboundary + lsign * l_coboundary
			else
				# add last (odd) edge to cycle
				rsign = rindex==EV[r_edge][1] ? 1 : -1
				cycle = cycle + rsign * r_coboundary
			end
			chain = setdiff(chain, findnz(cycle)[1])
		end
		for e in findnz(cscFE[f,:])[1]
			cscFE[f,e] = cycle[e]
		end
	end
	if exterior && size(V,1)==2
		# put matrix in form: first row outer cell; with opposite sign )
		V = convert(Array{Float64,2},transpose(V))
		EV = convert(ChainOp, SparseArrays.transpose(boundary_1(EV)))

		outer = Arrangement.get_external_cycle(V::Points, cscEV::ChainOp,
			cscFE::ChainOp)
		copFE = [ -cscFE[outer:outer,:];  cscFE[1:outer-1,:];  cscFE[outer+1:end,:] ]
		# induce coherent orientation of matrix rows (see examples/orient2d.jl)
		for k=1:size(copFE,2)
			spcolumn = findnz(copFE[:,k])
			if sum(spcolumn[2]) != 0
				row = spcolumn[1][2]
				sign = spcolumn[2][2]
				copFE[row,:] = -sign * copFE[row,:]
			end
		end
		return copFE
	else
		return cscFE
	end
end

"""
	u_boundary_2(FV::Cells, EV::Cells)::ChainOp

Return the unsigned `boundary_2` operator `C_2` -> `C_1`.
"""
u_boundary_2(EV, FV) = (u_coboundary_1(FV, EV))'

"""
	u_boundary_3(CV::Cells, FV::Cells)::ChainOp

Return the unsigned `boundary_3` operator `C_3` -> `C_2`.
"""
u_boundary_3(CV, FV) = (u_coboundary_2(CV, FV))'

"""
	u_coboundary_2( CV::Cells, FV::Cells[, convex=true::Bool] )::ChainOp

Unsigned 2-coboundary matrix `∂_2 : C_2 -> C_3` from 2-chain to 3-chain space.
Compute algebraically the *unsigned* coboundary matrix `∂_2` from
characteristic matrices of `CV` and `FV`. Currently usable *only* with complexes of *convex* cells.

#	Examples

## First example

(1) Compute the *boundary matrix* for a block of 3-cells of size ``[32,32,16]``;

(2) compute and show the *boundary* 2-cell array `boundary_2D_cells` by decodifying the (`mod 2`) result of multiplication of  the *boundary_3 matrix* `∂_2'`, transpose of *unsigned  coboundary_2* matrix  times the coordinate vector of the ``total`` 3-chain.

```julia
julia> using SparseArrays, Plasm
julia> V,(_,_,FV,CV) = cuboidGrid([32,32,16], true)
julia> ∂_2 = u_coboundary_2( CV, FV)
julia> coord_vect_of_all_3D_cells  = ones(size(∂_2,1),1)
julia> coord_vect_of_boundary_2D_cells = ∂_2' * coord_vect_of_all_3D_cells .% 2
julia> out = coord_vect_of_boundary_2D_cells
julia> boundary_2D_cells = [ FV[f] for f in findnz(sparse(out))[1] ]
julia> hpc = Plasm.lar2exploded_hpc(V, boundary_2D_cells)(1.,1.,1.)
julia> Plasm.view(hpc)
```
## Second example example

Using the boundary matrix of the `32 x 32 x 16` "image block" (better if stored on disk)
compute the boundary 2-complex of a random sub-image inside the block.

```julia
julia> coord_vect_of_segment = [x>0.25 ? 1 : 0  for x in rand(size(∂_2,1)) ]
julia> out = ∂_2' * coord_vect_of_segment .% 2
julia> boundary_2D_cells = [ FV[f] for f in findnz(sparse(out))[1] ]
julia> hpc = Plasm.lar2exploded_hpc(V, boundary_2D_cells)(1.1,1.1,1.1)
julia> Plasm.view(hpc)
```

"""
function u_coboundary_2( CV::Cells, FV::Cells, convex=true::Bool)::ChainOp
	cscCV = characteristicMatrix(CV)
	cscFV = characteristicMatrix(FV)
	temp = cscCV * cscFV'
	I,J,value = Int64[],Int64[],Int8[]
	for j=1:size(temp,2)
		nverts = length(FV[j])
		for i=1:size(temp,1)
			if temp[i,j] == nverts
				push!(I,i)
				push!(J,j)
				push!(value,1)
			end
		end
	end
	cscCF = SparseArrays.sparse(I,J,value)
	if !convex
		@assert "not yet implemented: TODO!"
	end
	return cscCF
end

"""
	chaincomplex( W::Points, EW::Cells )::Tuple{Array{Cells,1},Array{ChainOp,1}}

Chain 2-complex construction from basis of 1-cells.

From the minimal input, construct the whole
two-dimensional chain complex, i.e. the bases for linear spaces C_1 and
C_2 of 1-chains and  2-chains, and the signed coboundary operators from
C_0 to C_1 and from C_1 to C_2.

# Example

```julia
julia> W =
[0.0  0.0  0.0  0.0  1.0  1.0  1.0  1.0  2.0  2.0  2.0  2.0  3.0  3.0  3.0  3.0
0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0  0.0  1.0  2.0  3.0]
# output
2×16 Array{Float64,2}: ...

julia> EW =
[[1, 2],[2, 3],[3, 4],[5, 6],[6, 7],[7, 8],[9, 10],[10, 11],[11, 12],[13, 14],
[14, 15],[15, 16],[1, 5],[2, 6],[3, 7],[4, 8],[5, 9],[6, 10],[7, 11],[8, 12],
[9, 13],[10, 14],[11, 15],[12, 16]]
# output
24-element Array{Array{Int64,1},1}: ...

julia> V,bases,coboundaries = chaincomplex(W,EW)

julia> bases[1]	# edges
24-element Array{Array{Int64,1},1}: ...

julia> bases[2] # faces -- previously unknown !!
9-element Array{Array{Int64,1},1}: ...

julia> coboundaries[1] # coboundary_1
24×16 SparseMatrixCSC{Int8,Int64} with 48 stored entries: ...

julia> Matrix(coboundaries[2]) # coboundary_1: faces as oriented 1-cycles of edges
9×24 Array{Int8,2}:
-1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  0  0
0 -1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0  0
0  0 -1  0  0  1  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0  0  0
0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0  0
0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0  0
0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  1 -1  0  0  0  0
0  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1  0
0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1  0  0
0  0  0  0  0  0  0  0 -1  0  0  1  0  0  0  0  0  0  0  0  0  0  1 -1
```
"""
function chaincomplex( W, EW )
	V = convert(Array{Float64,2},LinearAlgebra.transpose(W))
	EV = convert(ChainOp, SparseArrays.transpose(boundary_1(EW)))

	V,cscEV,cscFE = planar_arrangement(V,EV)

	ne,nv = size(cscEV)
	nf = size(cscFE,1)
	EV = [findall(!iszero, cscEV[e,:]) for e=1:ne]
	FV = [collect(Set(vcat([EV[e] for e in findall(!iszero, cscFE[f,:])]...)))  for f=1:nf]

	function ord(cells)
		return [sort(cell) for cell in cells]
	end
	temp = copy(convert(ChainOp, LinearAlgebra.transpose(cscEV)))
	for k=1:size(temp,2)
		h = findall(!iszero, temp[:,k])[1]
		temp[h,k] = -1
	end
	cscEV = convert(ChainOp, LinearAlgebra.transpose(temp))
	bases, coboundaries = (ord(EV),ord(FV)), (cscEV,cscFE)
	return V',bases,coboundaries
end

"""
	chaincomplex( W::Points, FW::Cells, EW::Cells )
		::Tuple{ Array{Cells,1}, Array{ChainOp,1} }

Chain 3-complex construction from bases of 2- and 1-cells.

From the minimal input, construct the whole
two-dimensional chain complex, i.e. the bases for linear spaces C_1 and
C_2 of 1-chains and  2-chains, and the signed coboundary operators from
C_0 to C_1  and from C_1 to C_2.

# Example
```julia
julia> cube_1 = ([0 0 0 0 1 1 1 1; 0 0 1 1 0 0 1 1; 0 1 0 1 0 1 0 1],
[[1,2,3,4],[5,6,7,8],[1,2,5,6],[3,4,7,8],[1,3,5,7],[2,4,6,8]],
[[1,2],[3,4],[5,6],[7,8],[1,3],[2,4],[5,7],[6,8],[1,5],[2,6],[3,7],[4,8]] )

julia> cube_2 = Struct([t(0,0,0.5), r(0,0,pi/3), cube_1])

julia> V,FV,EV = struct2lar(Struct([ cube_1, cube_2 ]))

julia> V,bases,coboundaries = chaincomplex(V,FV,EV)

julia> (EV, FV, CV), (cscEV, cscFE, cscCF) = bases,coboundaries

julia> FV # bases[2]
18-element Array{Array{Int64,1},1}:
[1, 3, 4, 6]
[2, 3, 5, 6]
[7, 8, 9, 10]
[1, 2, 3, 7, 8]
[4, 6, 9, 10, 11, 12]
[5, 6, 11, 12]
[1, 4, 7, 9]
[2, 5, 11, 13]
[2, 8, 10, 11, 13]
[2, 3, 14, 15, 16]
[11, 12, 13, 17]
[11, 12, 13, 18, 19, 20]
[2, 3, 13, 17]
[2, 13, 14, 18]
[15, 16, 19, 20]
[3, 6, 12, 15, 19]
[3, 6, 12, 17]
[14, 16, 18, 20]

julia> CV # bases[3]
3-element Array{Array{Int64,1},1}:
[2, 3, 5, 6, 11, 12, 13, 14, 15, 16, 18, 19, 20]
[2, 3, 5, 6, 11, 12, 13, 17]
[1, 2, 3, 4, 6, 7, 8, 9, 10, 11, 12, 13, 17]

julia> cscEV # coboundaries[1]
34×20 SparseMatrixCSC{Int8,Int64} with 68 stored entries: ...

julia> cscFE # coboundaries[2]
18×34 SparseMatrixCSC{Int8,Int64} with 80 stored entries: ...

julia> cscCF # coboundaries[3]
4×18 SparseMatrixCSC{Int8,Int64} with 36 stored entries: ...
```
"""
function chaincomplex(V,FV,EV)
	W = convert(Points, copy(V)');
    cop_EV = coboundary_0(EV::Cells);
    cop_FE = coboundary_1(V, FV::Cells, EV::Cells);

    W, copEV, copFE, copCF = Arrangement.spatial_arrangement( W::Points, cop_EV::ChainOp, cop_FE::ChainOp)
	ne,nv = size(copEV)
	nf = size(copFE,1)
	nc = size(copCF,1)
	EV = [findall(!iszero, copEV[e,:]) for e=1:ne]
	FV = [collect(Set(vcat([EV[e] for e in findall(!iszero, copFE[f,:])]...)))  for f=1:nf]
	CV = [collect(Set(vcat([FV[f] for f in findall(!iszero, copCF[c,:])]...)))  for c=2:nc]
	function ord(cells)
		return [sort(cell) for cell in cells]
	end
	temp = copy(convert(ChainOp, LinearAlgebra.transpose(copEV)))
	for k=1:size(temp,2)
		h = findall(!iszero, temp[:,k])[1]
		temp[h,k] = -1
	end
	copEV = convert(ChainOp, LinearAlgebra.transpose(temp))
	bases, coboundaries = (ord(EV),ord(FV),ord(CV)), (copEV,copFE,copCF)
	W = convert(Points, (LinearAlgebra.transpose(V')))
	return W,bases,coboundaries
end

# Collect LAR models in a single LAR model
function collection2model(collection)
	W,FW,EW = collection[1]
	shiftV = size(W,2)
	for k=2:length(collection)
		V,FV,EV = collection[k]
		W = [W V]
		FW = [FW; FV + shiftV]
		EW = [EW; EV + shiftV]
		shiftV = size(W,2)
	end
	return W,FW,EW
end

function arrange2D(V,EV)
	cop_EV = coboundary_0(EV::Cells)
	cop_EW = convert(ChainOp, cop_EV)
	W = convert(Points,V')
	V, copEV, copFE = Arrangement.planar_arrangement(W::Points, cop_EW::ChainOp)
	EVs = FV2EVs(copEV, copFE) # polygonal face fragments

	triangulated_faces = triangulate2D(V, [copEV, copFE])
	FVs = convert(Array{Cells}, triangulated_faces)
	V = convert(Points,V')
	return V,FVs,EVs, copEV, copFE
end

"""
	pols2tria(W::Points, copEV::ChainOp,
			copFE::ChainOp, copCF::ChainOp)

Take a chain 3-complex and return arrays of boundary simplicial complexes.
Input and ouput Points (embedding geometry) are by columns.
Arrays of simplicial complexes are by body, by face, and by boundary of face.
```

```
"""
function pols2tria(W, copEV, copFE, copCF) # W by columns
	V = convert(Points,W')
	triangulated_faces = triangulate(V, [copEV, copFE])
	EVs = FV2EVs(copEV, copFE) # polygonal face fragments
	# triangulated_faces = [ item for (k,item) in enumerate(triangulated_faces)
	# 	if isdefined(triangulated_faces,k) && item ≠ Any[]  ]
	FVs = convert(Array{Cells}, triangulated_faces)
	CVs = []
	for cell in 1:copCF.m
		obj = []
        for f in copCF[cell, :].nzind
            triangles = triangulated_faces[f]
			append!(obj, triangles)
        end
		push!(CVs,obj)
    end
	V = convert(Points,V')
	return V,CVs,FVs,EVs
end
## Fs is the signed coord vector of a subassembly
## the logic is to compute the corresponding reduced coboundary matrices
## and finally call the standard method of the function.
function pols2tria(W, copEV, copFE, copCF, Fs) # W by columns
	# make copies of coboundary operators

	# compute the reduced copCF
	CFtriples = findnz(copCF)
	triples = [triple for triple in zip(CFtriples...)]
	newtriples = [(row,col,val) for (row,col,val) in triples if Fs[col] ≠ 0]
	newF = [k for (k,f) in enumerate(Fs) if Fs[k] ≠ 0]
	fdict = Dict( zip(newF, 1:length(newF)))
	triples = hcat([[row,fdict[col],val] for (row,col,val) in newtriples]...)
	newCF = sparse( triples[1,:], triples[2,:], triples[3,:] )
	copCF = convert( SparseMatrixCSC{Int8,Int64}, newCF )

	# compute the reduced copFE
	FEtriples = findnz(copFE)
	triples = [triple for triple in zip(FEtriples...)]
	newtriples = [(row,col,val) for (row,col,val) in triples if Fs[row] ≠ 0]
	newF = [k for (k,f) in enumerate(Fs) if Fs[k] ≠ 0]
	newcol = collect(Set([col for (row,col,val) in newtriples]))
	facedict = Dict( zip(newF, 1:length(newF)))
	edgedict = Dict( zip(newcol, 1:length(newcol)))
	triples = hcat([ [facedict[row],edgedict[col],val] for (row,col,val) in newtriples]...)
	newFE = sparse( triples[1,:], triples[2,:], triples[3,:] )
	copFE = convert( SparseMatrixCSC{Int8,Int64}, newFE )

	# compute the reduced copEV
	EVtriples = findnz(copEV)
	triples = [triple for triple in zip(EVtriples...)]
	newtriples = [(row,col,val) for (row,col,val) in triples if row in keys(edgedict)]
	# newcol = collect(Set([col for (row,col,val) in newtriples]))
	# vertdict = Dict( zip(newcol, 1:length(newcol)))
	# triples = hcat([[edgedict[row],vertdict[col],val] for (row,col,val) in newtriples]...)
	triples = hcat([[edgedict[row],col,val] for (row,col,val) in newtriples]...)
	newEV = sparse( triples[1,:], triples[2,:], triples[3,:] )
	copEV = convert( SparseMatrixCSC{Int8,Int64}, newEV )

	#W = convert(Points,W') # BOH...!!
	# finally compute the cells, faces, and edges of subassembly
	V,CVs,FVs,EVs = pols2tria(W, copEV, copFE, copCF)
	return V,CVs,FVs,EVs
end

"""
	permutationOrbits(perm::OrderedDict)::Array{Array{Int64,1},1}

Compute the ``cycles`` of a `perm` ``permutation`` of the first integers (starting from 1).

The `perm` parameter is an ``ordered dictionary``. The output is an
array of arrays of integers (``orbits``).

# Examples

```
julia> dict(List) = OrderedDict((i,x) for (i,x) in enumerate(List))
dict (generic function with 1 method)

julia> perm = dict([2, 3, 4, 5, 6, 7, 8, 1]);

julia> permutationOrbits(perm)
1-element Array{Array{Int64,1},1}:
 [2, 3, 4, 5, 6, 7, 8, 1, 2]

julia> perm = dict([3,9,8,12,10,7,2,11,6,4,1,5]);

julia> permutationOrbits(perm)
3-element Array{Array{Int64,1},1}:
 [3, 8, 11, 1, 3]
 [9, 6, 7, 2, 9]
 [12, 5, 10, 4, 12]

julia> permutationOrbits(Dict())
0-element Array{Array{Int64,1},1}

julia> permutationOrbits(Dict(1=>1))
1-element Array{Array{Int64,1},1}:
 [1, 1]
```
"""
function permutationOrbits(perm::OrderedDict)
	out = Array{Int64,1}[]
    while perm ≠ Dict()
        x = collect(keys(perm))[1]
        orbit = Int64[]
        while x in keys(perm)
            append!(orbit, perm[x])
			y,x = x,perm[x]
            delete!(perm,y)
		end
        append!(out, [ push!(orbit,orbit[1]) ]  )
	end
    return out
end

## Generation of oriented boundary faces
"""
	faces2polygons(copEV::ChainOp, copFE::ChainOp)::Array{Array{Array{Int64,1},1},1}

Generate an array of faces (array of polygons, given as oriented cycles of vertex indices).

The input is a chain of two ``chain operators`` (coboundaries).
The output is a list of list of polygons, given as oriented ``cycles of vertex`` indices.
``Outer`` polygons are counterclockwise-oriented; ``inner`` polygons are clockwise-oriented.
The first polygon of each face is the outer one.

#	Example
Both a polygon with ``holes`` and polygons ``fitting the holes`` are defined in the following.
The reader should note the ``mutual orientation`` of holes and polygons fitting them.
They are correctly derived from a sparse `copFE` matrix created by ``TGW algorithm`` in 2D.

```
julia> V = [540.313 540.313 2038.65 2038.65 1990.25 1990.25 583.951 583.951 2038.65 2038.65 1551.97 1990.25 1551.97 1990.25 409.035 2265.24 409.035 2265.24 540.313 540.313 583.951 583.951 934.346 1246.02 934.346 1246.02; -2129.59 -1653.96 -2129.59 -1493.96 -2064.15 -1493.96 -2064.15 -1653.96 -1104.9 -1221.14 -1104.9 -1137.88 -1137.88 -1221.14 -1007.31 -1007.31 -2210.65 -2210.65 -1109.81 -1360.09 -1104.1 -1360.09 -1104.9 -1104.9 -1137.88 -1137.88]

julia> copEV = sparse([1, 2, 1, 8, 2, 3, 3, 7, 4, 5, 4, 7, 5, 6, 6, 8, 9, 10, 9, 13, 10, 14, 11, 12, 11, 14, 12, 13, 15, 16, 15, 17, 16, 18, 17, 18, 19, 21, 19, 22, 20, 21, 20, 22, 23, 24, 23, 26, 24, 25, 25, 26], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16, 17, 17, 18, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 25, 25, 26, 26], Int8[-1, -1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, 1, -1, 1, -1, -1, -1, 1, 1, 1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, -1, 1, -1, -1, 1, 1, 1, -1, -1, 1, -1, 1, -1, 1, 1])

julia> copFE = sparse([1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 1, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 3, 3, 3, 3, 3, 4, 3, 4, 3, 4, 3, 4, 3, 5, 3, 5, 3, 5, 3, 5], [1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 16, 17, 18, 19, 19, 20, 20, 21, 21, 22, 22, 23, 23, 24, 24, 25, 25, 26, 26], Int8[-1, 1, 1, -1, 1, -1, -1, 1, 1, -1, 1, -1, 1, -1, -1, 1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, -1, 1, 1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 1, 1, -1])

julia> faces2polygons(copEV,copFE)
5-element Array{Array{Array{Int64,1},1},1}:
 [[1, 3, 4, 6, 5, 7, 8, 2]]
 [[9, 11, 13, 12, 14, 10]]
 [[2, 8, 7, 5, 6, 4, 3, 1], [10, 14, 12, 13, 11, 9], [15, 17, 18, 16], [19, 21, 22, 20], [24, 26, 25, 23]]
 [[20, 22, 21, 19]]
 [[23, 25, 26, 24]]

julia> model = (V, [ [[v] for v=1:size(V,2)], cop2lar(copEV) ])
([540.313 540.313 … 934.346 1246.02; -2129.59 -1653.96 … -1137.88 -1137.88], Array{Array{Int64,1},1}[[[1], [2], [3], [4], [5], [6], [7], [8], [9], [10]  …  [17], [18], [19], [20], [21], [22], [23], [24], [25], [26]], [[1, 2], [1, 3], [3, 4], [5, 6], [5, 7], [7, 8], [4, 6], [2, 8], [9, 10], [9, 11]  …  [16, 18], [17, 18], [19, 20], [21, 22], [19, 21], [20, 22], [23, 24], [23, 25], [25, 26], [24, 26]]])

julia> GL.VIEW( GL.numbering(200)( model, GL.COLORS[1], 0.1 ) );
```
"""
function faces2polygons(copEV,copFE)
	polygons = Array{Array{Int64,1},1}[]
	cycles = Array{Array{Array{Int64,1},1},1}[]
	for f=1:size(copFE,1)
		edges,signs = findnz(copFE[f,:])
		permutationMap = OrderedDict([ s>0 ? findnz(copEV[e,:])[1] : reverse(findnz(copEV[e,:])[1])
				for (e,s) in zip(edges,signs)])
		orbits = permutationOrbits(permutationMap)
		edgecycles = [[[ orbit[k], orbit[k+1] ] for k=1:length(orbit)-1]  for orbit in orbits]
		push!(polygons, [orbit[1:end-1] for orbit in orbits])
		push!(cycles, edgecycles)
	end
	return polygons,cycles
end


text:
# Autogenerated wrapper script for GLFW_jll for i686-w64-mingw32
export libglfw

using Libglvnd_jll
using Xorg_libXcursor_jll
using Xorg_libXi_jll
using Xorg_libXinerama_jll
using Xorg_libXrandr_jll
JLLWrappers.@generate_wrapper_header("GLFW")
JLLWrappers.@declare_library_product(libglfw, "glfw3.dll")
function __init__()
    JLLWrappers.@generate_init_header(Libglvnd_jll, Xorg_libXcursor_jll, Xorg_libXi_jll, Xorg_libXinerama_jll, Xorg_libXrandr_jll)
    JLLWrappers.@init_library_product(
        libglfw,
        "bin\\glfw3.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# charlm of lexicon. cluster all hidden states.

using TSne, Statistics

include("charlm.jl")

# load model
if !@isdefined(model)
    include(Knet.dir("src/models21/deepmap.jl"))
    include(Knet.dir("src/models21/fileio.jl"))
    model = charlm_model()
    charlm_loss(model, tst) # necessary to initialize weights, fix this
    setweights!(model, "charlm.jld2")
end
#@show charlm_loss(model,tst)

# load lexical forms for labels
if !@isdefined(lexical_forms)
    lexical_forms = last.(collect(trmor_data))
    lexical_vocab = trmor_data.tgt.vocab
    ltst = lexical_forms[rperm[1:80]]
    lval = lexical_forms[rperm[81:160]]
    ltrn = lexical_forms[rperm[161:end]]
end
#@show lexical_vocab.i2w[ltst[1][1,:]]
#@show surface_vocab.i2w[tst[1][1,:]]

# try clustering all hidden states
# model(charlm_input(tst[1])) |> "37×128×8 CUDA.CuArray{Float32, 3, CUDA.Mem.DeviceBuffer}"
# model[1:end-1](charlm_input(tst[1])) |> "350×128×8 CUDA.CuArray{Float32, 3, CUDA.Mem.DeviceBuffer}"

hidden_states(model, batch) = model[1:end-1](charlm_input(batch)) |> Array{Float64} # Warning: Float32 makes tsne much slower!
# hidden_states(model, batch{Int64}(B,T)) => (H,B,T+1)

function all_hidden_states(model, data::Vector{Matrix{Int64}})
    hs = map(data) do batch
        hbt = hidden_states(model, batch)
        htb = permutedims(hbt, (1,3,2)) # we need (H,T+1,B) before hcat to get hiddens for a word contiguous
        return reshape(htb, size(htb,1), :)
    end
    reduce(hcat, hs)
end

#@show summary(final_hidden_states(model, tst))

function insert_markers(word::String)
    # produce labels for all positions in the word
    chars = collect(word)
    map(0:length(chars)) do i
        String(vcat(chars[1:i], ['.'], chars[i+1:end]))
    end
end

function all_surface_labels(batch)
    words = reduce(*, surface_vocab.i2w[batch]; dims=2)
    reduce(vcat, insert_markers.(words))
end

all_surface_labels(data::Vector{Matrix{Int64}}) = reduce(vcat, all_surface_labels.(data))

if !@isdefined tstahs
    tstahs = all_hidden_states(model, tst) |> permutedims # (104452, 350)
    tstalb = all_surface_labels(tst)
end

if !@isdefined tst_ahs_tsne
    n = 3000
    r = randperm(size(tstahs,1))[1:n]
    tst_ahs_tsne = tsne(tstahs[r,:])
    tst_ahs_tsne_labels = tstalb[r]
end

using Plots
plotly()
scatter(tst_ahs_tsne[:,1], tst_ahs_tsne[:,2]; hover=tst_ahs_tsne_labels)

# From: https://github.com/nassarhuda/JuliaTutorials/blob/master/TSNE/TSNE.ipynb
# function my_tsne_plot(Y,d,labels; marker=5)
#     colorsvec = get(ColorSchemes.tab10, labels, :extrema)
#     plot(legend=false,axis=false,grid=false)
#     tsneplot = scatter!(Y[:,1], Y[:,2],
#                         marker=(marker,marker,:auto,stroke(0)),alpha=0.8,
#                         background_color=:black,
#                         color=colorsvec)
#     for i = 0:9
#         mx,my = mean(Y[findall(labels.==i),:],dims=1)
#         s = Shape(mx.+[-d,d,d,-d],my.+[-d,-d,d,d])
#         plot!(s,color=ColorSchemes.tab10[i+1],linewidth=0,alpha=0.9)
#         annotate!(mx,my,text(string(i),color=:white))
#     end
#     return tsneplot
# end

# function tsne_mnist(N=3000; plot=false)
#     # using MLDatasets, ColorSchemes, Plots
#     alldata, labels = MNIST.traindata(Float64)
#     data = reshape(permutedims(alldata[:, :, 1:N], (3, 1, 2)),
#                    N, size(alldata, 1)*size(alldata, 2))
#     labels = labels[1:N]
#     rescale(A; dims=1) = (A .- mean(A, dims=dims)) ./ max.(std(A, dims=dims), eps())
#     # Normalize the data, this should be done if there are large scale differences in the dataset
#     X = rescale(data, dims=1)
#     @show summary(X)
#     Y1 = tsne(X, 2, 50, 1000, 20.0)
#     plot ? my_tsne_plot(Y1, 5, labels) : Y1
# end

# Timing: tsne(tstfhs[1:1000,:], 2, 50, 1000, 20.0);

# mnist tsne:        charlm tsne: (why is it taking so much longer with Float32?)
# N    secs  KL      secs  KL
# 1000 3     0.7470  3     0.3175
# 2000 14    1.0122  13    0.3681
# 3000 31    1.1529  25    0.4125
# 4000 48    1.2079

nothing

max_stars_repo_name: denizyuret/research

max_issues_repo_name: denizyuret/research

max_forks_repo_name: denizyuret/research


text:
"""
    function tb_train!(mdl, opti, trn, vld=nothing; epochs=1,
                      lr_decay=nothing, lrd_freq=1, lrd_linear=false,
                      l2=0.0,
                      eval_size=0.2, eval_freq=1,
                      acc_fun=nothing,
                      mb_loss_freq=100,
                      cp_freq=nothing, cp_dir="checkpoints",
                      tb_dir="logs", tb_name="run",
                      tb_text=\"\"\"Description of tb_train!() run.\"\"\",
                      opti_args...)

Train function with TensorBoard integration. TB logs are written with
the TensorBoardLogger.jl package.
The model is updated (in-place) and the trained model is returned.

### Arguments:
+ `mdl`: model; i.e. forward-function for the net
+ `opti`: Knet-stype optimiser type
+ `trn`: training data; iterator to provide (x,y)-tuples with
        minibatches
+ `vld`: validation data; iterator to provide (x,y)-tuples with
        minibatches. Set to `nothing`, if not defined.

### Keyword arguments:
#### Optimiser:
+ `epochs=1`: number of epochs to train
+ `lr_decay=nothing`: do a leraning rate decay if not `nothing`:
        the value given is the final learning rate after `epochs*lrd_freq`
        steps of decay. `lr_decay` is only applied if both start learning rate
        `lr` and final learning rate `lr_decay` are defined explicitly.
        Example: `lr=0.01, lr_decay=0.001` will reduce the lr from
        0.01 to 0.001 during the training.
+ `lrd_freq=1`: frequency of learning rate decay steps. Default is
        to modify the lr after every epoch
+ `lrd_linear=false`: type of learning rate decay;
        If `false`, lr is modified
        by a constant factor (e.g. 0.9) resulting in an exponential decay.
        If `true`, lr is modified by the same step size, i.e. linearly.
+ `l2=0.0`: L2 regularisation; implemented as weight decay per
        parameter
+ `opti_args...`: optional keyword arguments for the optimiser can be specified
        (i.e. `lr`, `gamma`, ...).

#### Model evaluation:
+ `eval_size=0.2`: fraction of validation data to be used for calculating
        loss and accuracy for train and validation data during training.
+ `eval_freq=1`: frequency of evaluation; default=1 means evaluation is
        calculated after each epoch. With eval_freq=10 eveluation is
        calculated 10 times per epoch.
+ `acc_fun=nothing`: function to calculate accuracy. The function
        must implement the following signature: `fun(model; data)` where
        data is an iterator that provides (x,y)-tuples of minibatches.
        For classification tasks, `accuracy` from the Knet package is
        a good choice. For regression a correlation or mean error
        may be used.
+ `mb_loss_freq=100`: frequency of training loss reporting. default=100
        means that 100 loss-values per epoch will be logged to TensorBoard.
        If mb_loss_freq is greater then the number of minibatches,
        loss is logged for each minibatch.
+ `cp_freq=nothing`: frequency of model checkpoints written to disk.
        Default is `nothing`, i.e. no checkpoints are written.
        To write the model after each epoch with
        name `model` use freq=1; to write every 2 epochs freq=0.5.
+ `cp_dir="checkpoints"`: directory for checkpoints

#### TensorBoard:
TensorBoard log-directory is created from 3 parts:
`tb_dir/tb_name/<current date time>`.

+ `tb_dir="logs"`: root directory for tensorborad logs.
+ `tb_name="run"`: name of training run. `tb_name` will be used as
        directory name and should not include whitespace
+ `tb_text`:  description
        to be included in the TensorBoard log as *text* log.
"""
function tb_train!(mdl, opti, trn, vld=nothing; epochs=1,
                  lr_decay=nothing, lrd_freq=1, lrd_linear=false,
                  l2=0.0,
                  eval_size=0.2, eval_freq=1, acc_fun=nothing,
                  mb_loss_freq=100,
                  cp_freq=nothing, cp_dir="checkpoints",
                  tb_dir="logs", tb_name="run",
                  tb_text="""Description of tb_train!() run.""",
                  opti_args...)

    # use every n-th mb for evaluation (based on vld if defined):
    #
    n_trn = length(trn)

    if vld == nothing
        n_vld = 0
        eval_vld = nothing
        n_eval = Int(ceil(n_trn * eval_size))
    else
        n_vld = length(vld)
        n_eval = Int(ceil(n_vld * eval_size))
        nth_vld = cld(n_vld, n_eval)
        eval_vld = takenth(vld, nth_vld)
    end
    nth_trn = cld(n_trn, n_eval)
    eval_trn = takenth(trn, nth_trn)

    # frequencies of actions during training:
    #
    eval_nth = cld(n_trn, eval_freq)
    mb_loss_nth = cld(n_trn, mb_loss_freq)
    lr_nth = cld(n_trn, lrd_freq)
    if cp_freq != nothing
        cp_nth = cld(n_trn, cp_freq)
    end

    # do lr-decay only if lr is explicitly defined:
    #
    if lr_decay != nothing && haskey(opti_args, :lr)
       lr_decay = calc_d_η(opti_args[:lr], lr_decay, lrd_linear, lrd_freq*epochs)
    else
        lr_decay = nothing
    end

    # mk log directory:
    #
    start_time = Dates.now()
    tb_log_dir = joinpath(pwd(), tb_dir, tb_name,
                    Dates.format(start_time, "yyyy-mm-ddTHH-MM-SS"))
    println("Training $epochs epochs with $n_trn minibatches/epoch")
    if vld != nothing
        println("    (and $n_vld validation mbs).")
    end
    println("Evaluation is performed every $eval_nth minibatches (with $n_eval mbs).")
    println("Watch the progress with TensorBoard at: $tb_log_dir")

    # Tensorboard logger:
    #
    tbl = TensorBoardLogger.TBLogger(tb_log_dir,
                    min_level=Logging.Info)
    log_text(tbl, tb_log_dir, tb_name, start_time, tb_text,
             opti, trn, vld, epochs,
             lr_decay, lrd_freq, l2,
             cp_freq, opti_args)
    calc_and_report_loss(mdl, eval_trn, eval_vld, tbl, 0)

    if acc_fun != nothing
        calc_and_report_acc(mdl, acc_fun, eval_trn, eval_vld, tbl, 0)
    end

    # set optimiser - only if not yet set:
    #
    if first(params(mdl)).opt == nothing
        for p in params(mdl)
            p.opt = opti(;opti_args...)
        end
    else
        if haskey(opti_args, :lr)
            set_learning_rate(mdl, opti_args[:lr])
        end
    end

    # Training:
    #
    mb_losses = Float32[]
    @showprogress for (i, (x,y)) in enumerate(ncycle(trn,epochs))

        loss = @diff mdl(x,y)
        mb_loss = value(loss)

        if isnan(mb_loss)
            println("ERROR: training aborted because of loss value NaN!")
            break
        end

        for p in params(loss)
            Δw = grad(loss, p) + p .* l2
            # println("updating $i: $(p.opt.lr), Δw: -")
            Knet.update!(p, Δw)
        end

        # TensorBoard:
        #
        # println("mb_loss: $mb_loss"); flush(stdout)
        push!(mb_losses, mb_loss)
        if (i % eval_nth) == 0
            calc_and_report_loss(mdl, eval_trn, eval_vld, tbl, eval_nth)

            if acc_fun != nothing
                calc_and_report_acc(mdl, acc_fun, eval_trn, eval_vld,
                                    tbl, eval_nth)
            end
        end

        if (i % mb_loss_nth) == 0
            TensorBoardLogger.log_value(tbl,
                    "Minibatch loss (epoch = $n_trn steps)",
                    mean(mb_losses), step=i)
            # println("mb_loss-mean: $(mean(mb_losses))"); flush(stdout)
            mb_losses = Float32[]
        end

        # checkpoints:
        #
        if (cp_freq != nothing) && (i % cp_nth) == 0
            write_cp(mdl, i, tb_log_dir)
        end

        # lr decay:
        #
        if (lr_decay != nothing) && (i % lr_nth == 0)
            lr = first(params(mdl)).opt.lr
            lr = lrd_linear ? lr + lr_decay : lr * lr_decay
            @printf("Setting learning rate to η=%.2e\n", lr)
            set_learning_rate(mdl, lr)
            # for p in params(mdl)
            #     p.opt.lr = lr
            #     # println("adapting lr in $p to $(p.opt.lr)")
            # end
        end
    end

    println("Training finished with:")
    println("Training loss:       $(calc_loss(mdl, data=trn))")
    if acc_fun != nothing
        println("Training accuracy:   $(calc_acc(mdl, acc_fun, data=trn))")
    end

    if vld != nothing
        println("Validation loss:     $(calc_loss(mdl, data=vld))")
        if acc_fun != nothing
            println("Validation accuracy: $(calc_acc(mdl, acc_fun, data=vld))")
        end
    end

    # save final model:
    #
    if (cp_freq != nothing)
        write_cp(mdl, n_trn*epochs+1, tb_log_dir)
    end
    return mdl
end

function set_learning_rate(mdl, lr)

    for p in params(mdl)
        p.opt.lr = lr
    end
end

function log_text(tbl, tb_log_dir, tb_name, start_time, tb_text,
                  opti, trn, vld, epochs,
                  lr_decay, lrd_freq, l2,
                  cp_freq, opti_args)

    if vld == nothing
        vld = []
    end

    tb_log_text =
    "
NNHelferlein.jl tb_train!() log
" *
    "
 " *
    "
dir:  $tb_log_dir
 " *
    "
name: $tb_name
 " *
    "
time: $(Dates.format(start_time, "E, yyyy/mm/dd, HH:MM:SS"))
 " *
    "
 " *
    "
 " *
    "   $tb_text " *
    "
 " *
    "
Training parameters:
" *
    "
" *
    "
Training minibatches: $(length(trn))
" *
    "
Validation minibatches: $(length(vld))
" *
    "
Optimiser: $opti
" *
    "
Epochs: $epochs
"

    if lr_decay != nothing
        tb_log_text *= "
learning rate step size: $lr_decay with frequency $lrd_freq
"
    end
    if l2 > 0
        tb_log_text *= "
L2 regularisation: $l2
"
    end
    if cp_freq != nothing
        tb_log_text *= "
Checkpoints are saved $cp_freq times per epoch
"
    end

    for arg in keys(opti_args)
        tb_log_text *= "
$arg: $(opti_args[arg])
"
    end
    tb_log_text *= "
"

    TensorBoardLogger.log_text(tbl, "Description", tb_log_text, step=0)
    # with_logger(tbl) do
    #         @info "Description" text=tb_log_text log_step_increment=0
    # end
end

function write_cp(model, step, dir)

    dir_name = joinpath(dir, "checkpoints")
    if !isdir(dir_name)
        mkdir(dir_name)
    end
    fname = joinpath(dir_name, "checkpoint_$step.jld2")
    JLD2.@save fname model
end

# Helper to calc loss and acc with only ONE forward run:
#
function calc_loss(mdl; data)

    loss = 0.0
    for (x,y) in data
        loss += mdl(x,y)
    end
    return loss/length(data)
end

function calc_acc(mdl, fun; data)

    # acc = 0.0
    # for (x,y) in data
    #     acc += fun(mdl(x), y)
    # end
    # return acc/length(data)
    return fun(mdl, data=data)
end

# Helper for TensorBoardLogger:
#
function calc_and_report_loss(mdl, trn, vld, tbl, step)

    loss_trn = calc_loss(mdl, data=trn)

    if vld != nothing
        loss_vld = calc_loss(mdl, data=vld)
        with_logger(tbl) do
            @info "Evaluation Loss" train=loss_trn valid=loss_vld log_step_increment=step
        end
    else
        with_logger(tbl) do
            @info "Evaluation Loss" train=loss_trn log_step_increment=0
        end
    end
end

function calc_and_report_acc(mdl, acc_fun, trn, vld, tbl, step)

    acc_trn = calc_acc(mdl, acc_fun, data=trn)

    if vld != nothing
        acc_vld = calc_acc(mdl, acc_fun, data=vld)
        with_logger(tbl) do
            @info "Evaluation Accuracy" train=acc_trn valid=acc_vld log_step_increment=step
        end
    else
        with_logger(tbl) do
            @info "Evaluation Accuracy" train=acc_trn log_step_increment=step
        end
    end
end

# calc step size for lr-decay based on
# start rate, end rate and num steps:
#
function calc_d_η(η_start, η_end, lrd_linear, steps)
    if !lrd_linear
        d = log(η_end/η_start) / steps
        d = exp(d)
    else
        d = (η_end - η_start) / steps
    end
    return d
end

"""
    function predict_top5(mdl, x; top_n=5, classes=nothing)

Run the model `mdl` for data in `x` and print the top 5
predictions as softmax probabilities.

### Arguments:
+ `top_n`: print top *n* hits
+ `classes`: optional list of human readable class labels.
"""
function predict_top5(mdl, x; top_n=5, classes=nothing)

    y = predict(mdl, x, softmax=false)

    if classes == nothing
        classes = repeat(["-"], size(y)[1])
    end
    for (i,o) in enumerate(eachcol(y))
        o = Knet.softmax(vec(Array(o)))
        top = sortperm(vec(Array(o)), rev=true)[1:top_n]
        println("top-$top_n hits for sample $i: $top"); flush(stdout)

        @printf("%6s  %6s   %s\n", "softmax", "#", "class label")
        for t in top
            @printf(" %6.2f  %6i   \"%s\"\n", o[t], t, classes[t])
        end
        println(" ")
    end
    return y
end

"""
    function predict(mdl, x; softmax=false)

Return the prediction for x.

### Arguments:
+ `mdl`: executable network model
+ `x`: iterator providing minibatches
        of input data
+ `softmax`: if true and if model is a `::Classifier` the predicted
        softmax probabilities are returned instead of raw
        activations.
"""
function predict(mdl, x; softmax=false)

    if x isa AbstractArray
        y = mdl(x)
    else
        # ys = [mdl(i) for i in x]
        y = cat((mdl(i) for i in x)..., dims=2)
    end
    y = convert(Array{Float32}, y)

    if softmax && mdl isa Classifier
        return Knet.softmax(y)
    else
        return y
    end
end

"""
    function hamming_dist(p, t; accuracy=false, vocab=nothing, pad=0)

Return the Hamming distance between two sequences or two minibatches
of sequences. Predicted sequences `p` and teaching input sequences `t`
may be of different length but the number of sequences in the minibatch
must be the same.

### Arguments:
+ `p`, `t`: n-dimensional arrays of type `Int` with predictions
        and teaching input for a minibatch of sequences.
        Shape of the arrays must be identical except of the first dimension
        (i.e. the sequence length) that may differ between `p` and `t`.
+ `accuracy=false`: if `false`, the mean Hamming distance in the minibatch
        is returned (i.e. the average number of differences in the sequences).
        If `true`, the accuracy is returned
        for all not padded positions in a range (0.0 - 1.0).
+ `vocab=nothing`: target laguage vocabulary of type `NNHelferlein.WordTokenizer`.
        If defined,
        the padding token of `vocab` is used to mask all control tokens in the
        sequences (i.e. '<start>, <end>, <unknwon>, <pad>').
+ `pad=0`: if `vocab` is undefined, `pad` is used to pad `p`, if the sequence
        length of `p` is smaller than the length of `t`.
"""
function hamming_dist(p, t; accuracy=false, vocab=nothing, pad=0)

    # make 2d matrix of sequences:
    #
    n_seq_t = size(t)[1]
    t = reshape(t, n_seq_t,:)

    n_seq_p = size(p)[1]
    p = reshape(p, n_seq_p,:)

    n_mb = size(t)[2]

    # make all control-tokens the same:
    #
    if vocab != nothing
        START = vocab("<start>")
        END = vocab("<end>")
        UNK = vocab("<unknown>")
        PAD = vocab("<pad>")

        t[t .== START] .= PAD
        t[t .== END] .= PAD
        t[t .== UNK] .= PAD
    else
        PAD = pad
    end

    # make n_seq of p the same as t:
    #
    if n_seq_p > n_seq_t
        p = p[1:n_seq_t,:]
    end
    while size(p)[1] < n_seq_t
        p = vcat(p, repeat([PAD], inner=(1,n_mb)))
    end

    # mask preds same as teaching and count all
    # mask positions:
    #
    p[t .== PAD] .= PAD
    num_pad = length(t[t .== PAD])
    return accuracy ? (sum(p .== t) - num_pad)/(length(t)-num_pad) : sum(p .!= t)/n_mb
end


text:
module collision_ballWithBox_45Deg

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of SolidFileMesh
vmat2 = deepcopy(vmat1)                                           # material of convex decomposition of SolidFileMesh
vmat2.transparency = 0.7

cmat = "Steel"

@assembly ThreeDFiles begin
  world = Modia3D.Object3D(visualizeFrame=true)

  sphereMoving     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidSphere(0.5) , "Steel", vmat1; contactMaterial = cmat); r=[0.0, 0.0, 0.0],  fixed=false) # R=ModiaMath.rot2(-pi/2),

  box     = Modia3D.Object3D(world, Modia3D.Solid(Modia3D.SolidBox(1.0,4.0,10.0) , "Steel", vmat1; contactMaterial = cmat); r=[-3.0, 0.0, -2.0], fixed=true, R=ModiaMath.rot2(-pi/3))
end

gravField = Modia3D.UniformGravityField(g=9.81, n=[-1,0,0])
threeD = ThreeDFiles(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true, visualizeSupportPoints=true))

model = Modia3D.SimulationModel( threeD )
result = ModiaMath.simulate!(model; stopTime=1.3, tolerance=1e-8,interval=0.001, log=false)

println("... success of collision_ballWithBox_45Deg.jl!")

end


text:
APERTURES = [
    CircularAperture,
    CircularAnnulus,
    EllipticalAperture,
    EllipticalAnnulus,
    RectangularAperture,
    RectangularAnnulus
]

PARAMS = [
    (3),
    (3, 5),
    (3, 3, 0),
    (3, 5, 4, 0),
    (3, 5, 0),
    (3, 5, 4, 0)
]

###########################
# Some helpers for testing
area(ap::CircularAperture) = π * ap.r^2
area(ap::CircularAnnulus) = π * (ap.r_out^2 - ap.r_in^2)
area(ap::EllipticalAperture) = π * ap.a * ap.b
area(ap::EllipticalAnnulus) = π * ap.a_out * ap.b_out - π * ap.a_in * ap.b_in
area(ap::RectangularAperture) = ap.w * ap.h
area(ap::RectangularAnnulus) = ap.w_out * ap.h_out - ap.w_in * ap.h_in

@testset "outside - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = ones(10, 10)
    aperture = AP(-60, 60, params...)
    @test aperture_photometry(aperture, data).aperture_sum ≈ 0
end

@testset "inside zeros - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = zeros(40, 40)
    aperture = AP(20.0, 20.0, params...)

    table_cent = aperture_photometry(aperture, data, method = :center)
    table_sub = aperture_photometry(aperture, data, method = (:subpixel, 10))
    table_ex = aperture_photometry(aperture, data, method = :exact)

    @test table_ex.aperture_sum ≈ 0
    @test table_sub.aperture_sum ≈ 0
    @test table_cent.aperture_sum ≈ 0

end

@testset "inside ones - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = ones(40, 40)
    aperture = AP(20.0, 20.0, params...)

    table_cent = aperture_photometry(aperture, data, method = :center)
    table_sub = aperture_photometry(aperture, data, method = (:subpixel, 10))
    table_ex = aperture_photometry(aperture, data, method = :exact)

    true_flux = area(aperture)

    @test table_ex.aperture_sum ≈ true_flux
    @test table_sub.aperture_sum ≈ table_ex.aperture_sum atol = 0.1
    @test table_cent.aperture_sum ≤ table_ex.aperture_sum

end

@testset "interface" begin
    data = zeros(40, 40)
    err = zeros(40, 40)
    aperture = CircularAperture(20.0, 20.0, 5.0)

    t1 = aperture_photometry(aperture, data)
    t2 = aperture_photometry(aperture, data, err)

    # 1.0 compat (no hasproperty function)
    hasfunc = VERSION < v"1.1" ? haskey : hasproperty

    @test !hasfunc(t1, :aperture_sum_err)
    @test t2.aperture_sum_err == 0

    apertures = CircularAperture.(20, 20, [1, 2, 3])
    t1 = aperture_photometry(apertures, data)
    t2 = aperture_photometry(apertures, data, err)

    @test !hasfunc(t1, :aperture_sum_err)
    @test t2.aperture_sum_err == zeros(3)
end

@testset "type stability - $AP" for (AP, params) in zip(APERTURES, PARAMS)
    data = zeros(40, 40)
    err = zeros(40, 40)
    aperture = AP(20.0, 20.0, params...)

    @inferred aperture_photometry(aperture, data, method = :center)
    @inferred aperture_photometry(aperture, data, method = (:subpixel, 10))
    @inferred aperture_photometry(aperture, data, method = :exact)

    @inferred aperture_photometry(aperture, data, err, method = :center)
    @inferred aperture_photometry(aperture, data, err, method = (:subpixel, 10))
    @inferred aperture_photometry(aperture, data, err, method = :exact)
end

@testset "photometry - circular" begin
    function test_aperture(data, aperture)
        error = ones(size(data))

        table_cent = aperture_photometry(aperture, data, error, method = :center)
        table_sub = aperture_photometry(aperture, data, error, method = (:subpixel, 12))
        table_ex = aperture_photometry(aperture, data, error, method = :exact)

        true_flux = area(aperture)
        true_err = sqrt(true_flux)

        @test table_ex.aperture_sum ≈ true_flux
        @test table_sub.aperture_sum ≈ table_ex.aperture_sum atol = 0.1
        @test table_cent.aperture_sum < table_ex.aperture_sum

        @test table_ex.aperture_sum_err ≈ true_err
        @test table_sub.aperture_sum_err ≈ table_ex.aperture_sum_err atol = 0.1
        @test table_cent.aperture_sum_err < table_ex.aperture_sum_err
    end

    @testset "errors - CircularAperture" begin
        data = ones(40, 40)
        aperture = CircularAperture(20, 20, 10)
        test_aperture(data, aperture)
    end

    @testset "errors - CircularAnnulus" begin
        data = ones(40, 40)
        aperture = CircularAnnulus(20, 20, 8, 10)
        test_aperture(data, aperture)
    end

    @testset "partial overlap" begin
        data = ones(20, 20)
        error = ones(size(data))
        positions = [10.5 10.5; 1 1; 1 20; 20 1; 20 20]
        apertures = [CircularAperture(positions[i, :], 5) for i in axes(positions, 1)]

        table = aperture_photometry(apertures, data, error)
        @test table.aperture_sum[1] ≈ 25π
        @test all(table.aperture_sum[2:end] .< 25π)
    end
end # photometry - circular

@testset "photometry - elliptical" begin
    function test_elliptical_aperture(data, aperture)
        error = ones(size(data))

        table_ex = aperture_photometry(aperture, data, error, method = :exact)
        table_cent = aperture_photometry(aperture, data, error, method = :center)
        table_sub = aperture_photometry(aperture, data, error, method = (:subpixel, 128))

        true_flux = area(aperture)
        true_err = sqrt(true_flux)

        @test table_ex.aperture_sum ≈ true_flux
        @test table_sub.aperture_sum ≈ true_flux rtol = 1e-3
        @test table_cent.aperture_sum <= table_sub.aperture_sum

        @test table_ex.aperture_sum_err ≈ true_err
        @test table_sub.aperture_sum_err ≈ true_err rtol = 1e-3
        @test table_cent.aperture_sum_err <= true_err
    end

    @testset "errors - EllipticalAperture" begin
        data = ones(40, 40)
        aperture = EllipticalAperture(20, 20, 10, 10, 0)
        test_elliptical_aperture(data, aperture)

    end

    @testset "errors - EllipticalAnnulus" begin
        data = ones(40, 40)
        aperture = EllipticalAnnulus(20, 20, 8, 10, 10, 0)
        test_elliptical_aperture(data, aperture)
    end

    @testset "partial overlap elliptical aperture" begin
        data = ones(20, 20)
        error = ones(size(data))
        positions = [10.5 10.5; 1 1; 1 20; 20 1; 20 20]
        apertures = [EllipticalAperture(positions[i, :], 5, 5, 0) for i in axes(positions, 1)]

        table = aperture_photometry(apertures, data, error, method = (:subpixel, 128))
        @test table.aperture_sum[1] ≈ 25π rtol = 1e-3
        @test all(table.aperture_sum[2:end] .< 25π)
    end

end # photometry elliptical

@testset "photometry - rectangular" begin
    function test_aperture(data, aperture)
        error = ones(size(data))

        table_cent = aperture_photometry(aperture, data, error, method = :center)
        table_sub = aperture_photometry(aperture, data, error, method = (:subpixel, 64))

        true_flux = area(aperture)
        true_err = sqrt(true_flux)

        @test table_sub.aperture_sum ≈ true_flux rtol = 1e-2
        @test table_cent.aperture_sum < table_sub.aperture_sum

        @test table_sub.aperture_sum_err ≈ true_err rtol = 1e-2
        @test table_cent.aperture_sum_err < table_sub.aperture_sum_err
    end

    @testset "errors - RectangularAperture" begin
        data = ones(40, 40)
        aperture = RectangularAperture(20, 20, 10, 5, 0)
        test_aperture(data, aperture)
    end

    @testset "errors - RectangularAnnulus" begin
        data = ones(40, 40)
        aperture = RectangularAnnulus(20, 20, 8, 10, 4, 0)
        test_aperture(data, aperture)
    end

    @testset "partial overlap" begin
        data = ones(20, 20)
        error = ones(size(data))
        positions = [10.5 10.5; 1 1; 1 20; 20 1; 20 20]
        apertures = [RectangularAperture(positions[i, :], 10, 10, 0) for i in axes(positions, 1)]

        table = aperture_photometry(apertures, data, error, method = (:subpixel, 64))
        @test table.aperture_sum[1] ≈ 100 rtol = 1e-2
        @test all(table.aperture_sum[2:end] .< 100)
    end
end # photometry - circular


text:
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"

"""
Open all libraries
"""
function __init__()
    global ATK_artifact_dir = abspath(artifact"ATK")
    global Bzip2_artifact_dir = abspath(artifact"Bzip2")
    global Cairo_artifact_dir = abspath(artifact"Cairo")
    global Dbus_artifact_dir = abspath(artifact"Dbus")
    global Expat_artifact_dir = abspath(artifact"Expat")
    global Fontconfig_artifact_dir = abspath(artifact"Fontconfig")
    global FreeType2_artifact_dir = abspath(artifact"FreeType2")
    global FriBidi_artifact_dir = abspath(artifact"FriBidi")
    global GTK3_artifact_dir = abspath(artifact"GTK3")
    global Gettext_artifact_dir = abspath(artifact"Gettext")
    global Glib_artifact_dir = abspath(artifact"Glib")
    global Graphene_artifact_dir = abspath(artifact"Graphene")
    global Graphite2_artifact_dir = abspath(artifact"Graphite2")
    global HarfBuzz_artifact_dir = abspath(artifact"HarfBuzz")
    global JpegTurbo_artifact_dir = abspath(artifact"JpegTurbo")
    global LZO_artifact_dir = abspath(artifact"LZO")
    global Libepoxy_artifact_dir = abspath(artifact"Libepoxy")
    global Libffi_artifact_dir = abspath(artifact"Libffi")
    global Libgcrypt_artifact_dir = abspath(artifact"Libgcrypt")
    global Libglvnd_artifact_dir = abspath(artifact"Libglvnd")
    global Libgpg_error_artifact_dir = abspath(artifact"Libgpg_error")
    global Libiconv_artifact_dir = abspath(artifact"Libiconv")
    global Libtiff_artifact_dir = abspath(artifact"Libtiff")
    global Libuuid_artifact_dir = abspath(artifact"Libuuid")
    global PCRE_artifact_dir = abspath(artifact"PCRE")
    global Pango_artifact_dir = abspath(artifact"Pango")
    global Pixman_artifact_dir = abspath(artifact"Pixman")
    global Wayland_artifact_dir = abspath(artifact"Wayland")
    global Wayland_protocols_artifact_dir = abspath(artifact"Wayland_protocols")
    global X11_artifact_dir = abspath(artifact"X11")
    global XML2_artifact_dir = abspath(artifact"XML2")
    global XSLT_artifact_dir = abspath(artifact"XSLT")
    global Xorg_compositeproto_artifact_dir = abspath(artifact"Xorg_compositeproto")
    global Xorg_damageproto_artifact_dir = abspath(artifact"Xorg_damageproto")
    global Xorg_fixesproto_artifact_dir = abspath(artifact"Xorg_fixesproto")
    global Xorg_glproto_artifact_dir = abspath(artifact"Xorg_glproto")
    global Xorg_inputproto_artifact_dir = abspath(artifact"Xorg_inputproto")
    global Xorg_kbproto_artifact_dir = abspath(artifact"Xorg_kbproto")
    global Xorg_libX11_artifact_dir = abspath(artifact"Xorg_libX11")
    global Xorg_libXau_artifact_dir = abspath(artifact"Xorg_libXau")
    global Xorg_libXcomposite_artifact_dir = abspath(artifact"Xorg_libXcomposite")
    global Xorg_libXcursor_artifact_dir = abspath(artifact"Xorg_libXcursor")
    global Xorg_libXdamage_artifact_dir = abspath(artifact"Xorg_libXdamage")
    global Xorg_libXdmcp_artifact_dir = abspath(artifact"Xorg_libXdmcp")
    global Xorg_libXext_artifact_dir = abspath(artifact"Xorg_libXext")
    global Xorg_libXfixes_artifact_dir = abspath(artifact"Xorg_libXfixes")
    global Xorg_libXi_artifact_dir = abspath(artifact"Xorg_libXi")
    global Xorg_libXinerama_artifact_dir = abspath(artifact"Xorg_libXinerama")
    global Xorg_libXrandr_artifact_dir = abspath(artifact"Xorg_libXrandr")
    global Xorg_libXrender_artifact_dir = abspath(artifact"Xorg_libXrender")
    global Xorg_libXtst_artifact_dir = abspath(artifact"Xorg_libXtst")
    global Xorg_libpthread_stubs_artifact_dir = abspath(artifact"Xorg_libpthread_stubs")
    global Xorg_libxcb_artifact_dir = abspath(artifact"Xorg_libxcb")
    global Xorg_libxkbfile_artifact_dir = abspath(artifact"Xorg_libxkbfile")
    global Xorg_randrproto_artifact_dir = abspath(artifact"Xorg_randrproto")
    global Xorg_recordproto_artifact_dir = abspath(artifact"Xorg_recordproto")
    global Xorg_renderproto_artifact_dir = abspath(artifact"Xorg_renderproto")
    global Xorg_util_macros_artifact_dir = abspath(artifact"Xorg_util_macros")
    global Xorg_xcb_proto_artifact_dir = abspath(artifact"Xorg_xcb_proto")
    global Xorg_xextproto_artifact_dir = abspath(artifact"Xorg_xextproto")
    global Xorg_xineramaproto_artifact_dir = abspath(artifact"Xorg_xineramaproto")
    global Xorg_xkbcomp_artifact_dir = abspath(artifact"Xorg_xkbcomp")
    global Xorg_xkeyboard_config_artifact_dir = abspath(artifact"Xorg_xkeyboard_config")
    global Xorg_xproto_artifact_dir = abspath(artifact"Xorg_xproto")
    global Xorg_xtrans_artifact_dir = abspath(artifact"Xorg_xtrans")
    global Zlib_artifact_dir = abspath(artifact"Zlib")
    global Zstd_artifact_dir = abspath(artifact"Zstd")
    global adwaita_icon_theme_artifact_dir = abspath(artifact"adwaita_icon_theme")
    global at_spi2_atk_artifact_dir = abspath(artifact"at_spi2_atk")
    global at_spi2_core_artifact_dir = abspath(artifact"at_spi2_core")
    global gdk_pixbuf_artifact_dir = abspath(artifact"gdk_pixbuf")
    global hicolor_icon_theme_artifact_dir = abspath(artifact"hicolor_icon_theme")
    global iso_codes_artifact_dir = abspath(artifact"iso_codes")
    global libpng_artifact_dir = abspath(artifact"libpng")
    global xkbcommon_artifact_dir = abspath(artifact"xkbcommon")
end  # __init__()


text:
module TestDictArray

using FactCheck
using DataCubes
using DataCubes.setna!

facts("DictArray tests") do
  context("constructor tests") do
    context("basic tests") do
      col1 = nalift(rand(10, 50))
      col2 = nalift(reshape(1:500, 10, 50))
      col3 = nalift(reshape(map(i->string("sym_",i), 1:500), 10, 50))
      col4 = nalift(hcat(rand(10,30), fill(:testsym, 10, 20)))
      d1 = DictArray(c1=col1, c2=col2, c3=col3, c4=col4)
      d2 = DictArray(LDict([:c3, :c4, :c1, :c2], Any[col3, col4, col1, col2]), [:c1, :c2, :c3, :c4])
      d3 = DictArray(:c1=>col1, :c2=>col2, :c3=>col3, :c4=>col4)
      d4 = DictArray((:c1,col1), (:c2,col2), (:c3,col3), (:c4,col4))
      @fact d1 --> d2
      @fact d2 --> d3
      @fact d3 --> d4
      @fact darr(Any[LDict(:a=>Nullable(1),:b=>Nullable{Int}()),LDict(:a=>Nullable(3),:b=>Nullable(4))]) --> @darr(a=[1,3], b=[NA,4])
      @fact convert(DictArray, [LDict(:a=>Nullable(3),:b=>Nullable{Int}()), LDict(:a=>Nullable(5),:b=>Nullable(3))]) --> @darr(a=[3,5],b=[NA,3])
      @fact convert(DictArray, nalift([LDict(:a=>Nullable(3),:b=>Nullable{Int}()), LDict(:a=>Nullable(5),:b=>Nullable(3))])) --> @darr(a=[3,5],b=[NA,3])
      @fact convert(DictArray, @darr(1=>[1 2 3;4 5 6])) --> darr(1=>[1 2 3;4 5 6])
      #@fact_throws @darr([1,2,3])
      #@fact_throws @darr(Dict(:a=>3))
      #@fact_throws @darr(darr(a=[1,2,3]),darr(a=[1,2,3]))
      #@fact_throws @darr()
    end
    context("constructing DictArrays with null elements...") do
      col1 = nalift(rand(10, 50))
      setna!(col1, 5:10,1:30)
      col2 = nalift(reshape(1:500, 10, 50))
      setna!(col2, 1:7,40:45)
      col3 = nalift(reshape(map(i->string("str_",i), 1:500), 10, 50))
      setna!(col3, 2:10,1:20)
      col4 = nalift(hcat(rand(10,30), fill(:testsym, 10, 20)))
      setna!(col4, 5:10,30:40)
      d1 = DictArray(c1=col1, c2=col2, c3=col3, c4=col4)
      d2 = DictArray(LDict([:c3, :c4, :c1, :c2], Any[col3, col4, col1, col2]), [:c1, :c2, :c3, :c4])
      d3 = DictArray(:c1=>col1, :c2=>col2, :c3=>col3, :c4=>col4)
      d4 = DictArray((:c1,col1), (:c2,col2), (:c3,col3), (:c4,col4))
      @fact d1 --> d2
      @fact d2 --> d3
      @fact d3 --> d4
      @fact darr(darr(a=[1 2;3 4;5 6],b=["abc" 'a';1 2;:m "xyz"],:c=>[1.0 1.5;:sym 'a';"X" "Y"]), c=[1 2;3 4;5 6], :d=>map(Nullable, [1 2;3 4;5 6])) --> @darr(:c=>[1 2;3 4;5 6],:a=>[1 2;3 4;5 6],:b=>["abc" 'a';1 2;:m "xyz"],:d=>[1 2;3 4;5 6])
      @fact @darr(a=[1,2,3],b=1) --> darr(a=[1,2,3],b=[1,1,1])
    end
  end
  context("array related method tests") do
    d = DictArray(:a=>nalift(reshape(1:20, 5, 4)),
                  :b=>2.0*nalift(reshape(1:20, 5, 4)),
                  2=>nalift(fill(:sym, 5, 4)),
                  :third=>nalift(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    @fact d[3,2] --> LDict(:a=>Nullable(8), :b=>Nullable(16.0), 2=>Nullable(:sym), :third=>Nullable("str_8"))
    @fact d[1:3,2] --> DictArray(:a=>nalift([6,7,8]),:b=>nalift([12.0,14.0,16.0]),2=>nalift([:sym,:sym,:sym]),:third=>nalift(["str_6","str_7","str_8"]))
    @fact d[2:2,2:4] --> DictArray(:a=>nalift([7 12 17]),:b=>nalift([14.0 24.0 34.0]),2=>nalift([:sym :sym :sym]),:third=>nalift(["str_7" "str_12" "str_17"]))
    @fact findfirst(d, LDict(:a=>Nullable(8), :b=>Nullable(16.0), 2=>Nullable(:sym), :third=>Nullable("str_8"))) --> 8
    @fact (setindex!(d,LDict(:a=>Nullable(108),:b=>Nullable(216.0),2=>Nullable(:newsym),:third=>Nullable("newstr_8")), 3, 2);d[3,2]) -->
      LDict(:a=>Nullable(108),:b=>Nullable(216.0),2=>Nullable(:newsym),:third=>Nullable("newstr_8"))
    @fact (setindex!(d,Dict(:a=>Nullable(105),:b=>Nullable(210.0),2=>Nullable(:newsym2),:third=>Nullable("newstr_8")), 3, 2);d[3,2]) -->
      LDict(:a=>Nullable(105),:b=>Nullable(210.0),2=>Nullable(:newsym2),:third=>Nullable("newstr_8"))
    @fact (setindex!(d,(Nullable(1008),Nullable(2016.0),Nullable(:newsym3),Nullable("newstr_80")), 3, 2);d[3,2]) -->
      LDict(:a=>Nullable(1008),:b=>Nullable(2016.0),2=>Nullable(:newsym3),:third=>Nullable("newstr_80"))
    @fact (setindex!(d,
                     DictArray(:a=>nalift(reshape(collect(100:108),3,3)),
                     :b=>nalift(2.0*reshape(collect(100:108),3,3)),
                     2=>nalift(fill(:newval,3,3)),
                     :third=>nalift(fill("newstr",3,3))),
                     2:4, 1:3);d[2,1]) -->
      LDict(:a=>Nullable(100),:b=>Nullable(200.0),2=>Nullable(:newval),:third=>Nullable("newstr"))
    @fact size(d) --> (5,4)
    #@fact eltype(d) --> LDict{Any}
    #@fact eltype(pick(d,[:a,:third])) --> LDict{Symbol}
    @fact pick(d,:a)[1,1].value --> 1
    @fact pick(d,:b)[1,1].value --> 2.0
    @fact pick(d,(:a,))[1][1,1].value --> 1
    @fact pick(d,(:b,))[1][1,1].value --> 2.0
    @fact peel(d)[:a][2,2].value --> 103
    @fact peel(d)[:b][2,2].value --> 206.0
    @fact pick(d, 2) --> peel(d)[2]
    @fact endof(d) --> length(d)
    @fact transpose(d) --> DictArray(mapvalues(transpose, d.data))
    @fact permutedims(d, [1,2]) --> DictArray(mapvalues(v->permutedims(v,[1,2]), d.data))
    @fact DictArray(a=nalift([1,2,3]),b=nalift([1.0,2.0,3.0])) == DictArray(a=nalift([1,2,3]),b=nalift([1.0,2.0,3.0])) --> true
    @fact copy(d) !== d --> true
    @fact size(cat(1, d, d)) --> ntuple(n->n==1 ? size(d,1)*2 : size(d,n), ndims(d))
    @fact size(cat(2, d, d)) --> ntuple(n->n==2 ? size(d,2)*2 : size(d,n), ndims(d))
    @fact size(vcat(d, d)) --> ntuple(n->n==1 ? size(d,1)*2 : size(d,n), ndims(d))
    @fact size(hcat(d, d)) --> ntuple(n->n==2 ? size(d,2)*2 : size(d,n), ndims(d))
    @fact size(repeat(d, inner=[5,2], outer=[3,4])) --> (size(d,1)*5*3, size(d,2)*2*4)
    @fact @darr(a=[1 2 3;4 5 6], 'x'=[:a :b :c;:x :y :z])[:a] --> nalift([1 2 3;4 5 6])
    @fact @darr(a=1.0*[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])[:a,:b] --> Any[1.0*nalift([1 2 3;4 5 6]), nalift([:a :b :c;:x :y :z])]
    @fact @darr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])[[:a,:b]] --> @darr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])
    @fact @darr(a=[1 2 3;4 5 6], b=[:a :b :c;:x :y :z])[(:a,:b)] --> (nalift([1 2 3;4 5 6]), nalift([:a :b :c;:x :y :z]))
    @fact (t=darr(a=[1 2 3;4 5 6],b=[:a :b :c;:x :y :z]);t[1:2,1]=darr(a=[11,12],b=[:u,:v]);t) --> darr(a=[11 2 3;12 5 6],b=[:u :b :c;:v :y :z])
    @fact (t=darr(a=[1 2 3;4 5 6],b=[:a :b :c;:x :y :z]);DataCubes.setindex_nocheck!(t,darr(a=[11,12],b=[:u,:v]),1:2,1);t) --> darr(a=[11 2 3;12 5 6],b=[:u :b :c;:v :y :z])
    @fact (t=darr(a=[1 2 3;4 5 6],b=[:a :b :c;:x :y :z]);DataCubes.setindex_nocheck!(t,nalift(LDict(:a=>11,:b=>:u)),1,1);t) --> darr(a=[11 2 3;4 5 6],b=[:u :b :c;:x :y :z])
    @fact findfirst(darr(a=[1,2,3],b=['a','b','c']),(Nullable(2),Nullable('b'))) --> 2
  end
  context("additional method tests") do
    d1 = darr(:a=>nalift(reshape(1:20, 5, 4)),
              2=>nalift(fill(:sym, 5, 4)),
              :third=>nalift(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    d2 = darr(:x=>nalift(reshape(1:20, 5, 4)),
              :third=>nalift(fill(:sym, 5, 4)),
              :z=>nalift(reshape(map(i->string("str_",i), 1:20), 5, 4)))
    @fact pick(merge(d1, d2), [:a,2]) --> pick(d1, [:a,2])
    @fact pick(merge(d1, d2), [:x,:third,:z]) --> pick(d2, [:x,:third,:z])
    @fact merge(d1, d2).data.keys --> Any[:a,2,:third,:x,:z]
    @fact delete(merge(d1, d2), :a,2,:third) --> pick(d2, [:x, :z])
    @fact merge(darr(a=[1 2 3;4 5 6]), b=3.0) --> darr(a=[1 2 3;4 5 6], b=fill(3.0,2,3))
    @fact pick(d1, [:a]) --> DictArray(a=pick(d1, :a))
    @fact mapslices(x->x,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[]) --> @darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"])
    #@fact mapslices(x->x,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[1]) --> Any[@darr(a=[1,4], b=["a","d"]),@darr(a=[2,5], b=["b","e"]),@darr(a=[3,6],b=["c","f"])]
    @fact mapslices(x->x,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[1]) --> @darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"])
    @fact mapslices(x->x,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),1) --> @darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"])
    @fact mapslices(x->x[1],@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[1]) --> @darr(a=[1,2,3], b=["a","b","c"])
    @fact mapslices(x->x[1],@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[2]) --> @darr(a=[1,4], b=["a","d"])
    @fact mapslices(length,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[1]) --> @rap @nalift [2,2,2]
    @fact mapslices(length,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),[2]) --> @rap @nalift [3,3]
    @fact mapslices(length,@darr(a=[1 2 3;4 5 6],b=["a" "b" "c";"d" "e" "f"]),2) --> @rap @nalift [3,3]
    @fact mapslices(x->sum(map(x->x.value,peel(x)[:a])),@darr(a=[1,2,3]),[1]).value --> 6
    @fact mapslices(x->x,@larr(a=[1,2,3]),[1]) --> @larr(a=[1,2,3])
    @fact mapslices(x->x,@larr(a=[1,2,3]), [1]) --> mapslices(x->x,@larr(a=[1,2,3]), [])
    @fact mapslices(identity, @darr(a=[1 2 3;4 5 6]), []) --> @darr(a=[1 2 3;4 5 6])
    @fact mapslices(x->darr(b=x[:a]+10),darr(a=[10,9]),[1]) --> darr(b=[20,19])
    @fact mapslices(x->(x=x[:a].value;[x,x+1,x+10]),darr(a=[10,9]),[]) --> nalift([10 9;11 10;20 19])
    @fact mapslices(x->msum(x), darr(a=reshape(1:15,5,3), b=1.0*reshape(1:15,5,3)), [1]) --> @update(darr(a=[1 6 11;3 13 23;6 21 36;10 30 50;15 40 65]), b=1.0*_a)
    @fact mapslices(x->msum(x), darr(a=reshape(1:15,5,3), b=1.0*reshape(1:15,5,3)), [2]) --> update(darr(a=[1 2 3 4 5;7 9 11 13 15;18 21 24 27 30]), :b=>d->1.0*d[:a])

    @fact mapslices(x->darr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O]), darr(a=rand(3,5,2)), [1]) --> darr(l=reshape(repeat(collect(1:4),outer=[10]),4,5,2),p=reshape(repeat([:U,:V,:W,:O],outer=[10]),4,5,2))
    @fact mapslices(x->darr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O]), darr(a=rand(3,5,2)), [1]) --> darr(l=reshape(repeat(collect(1:4),outer=[10]),4,5,2),p=reshape(repeat([:U,:V,:W,:O],outer=[10]),4,5,2))
    @fact mapslices(x->larr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O]), darr(a=rand(3,5,2)), [1]) --> larr(l=reshape(repeat(collect(1:4),outer=[10]),4,5,2),p=reshape(repeat([:U,:V,:W,:O],outer=[10]),4,5,2))
    @fact mapslices(x->larr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O], axis1=[:m,:n,:p,:q]), darr(a=rand(3,5,2)), [1]) --> larr(l=reshape(repeat(collect(1:4),outer=[10]),4,5,2),p=reshape(repeat([:U,:V,:W,:O],outer=[10]),4,5,2),axis1=[:m,:n,:p,:q])
    @fact mapslices(x->[11,21,31,41], darr(a=rand(3,5,2)), [1]) --> reshape(repeat(nalift([11,21,31,41]),outer=[10]),4,5,2)

    @fact mapslices(x->darr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O]), darr(a=rand(3,5,2)), [2]) --> darr(l=reshape(repeat(collect(1:4),outer=[6]),4,3,2),p=reshape(repeat([:U,:V,:W,:O],outer=[6]),4,3,2))
    @fact mapslices(x->darr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O]), darr(a=rand(3,5,2)), [2]) --> darr(l=reshape(repeat(collect(1:4),outer=[6]),4,3,2),p=reshape(repeat([:U,:V,:W,:O],outer=[6]),4,3,2))
    @fact mapslices(x->larr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O]), darr(a=rand(3,5,2)), [2]) --> larr(l=reshape(repeat(collect(1:4),outer=[6]),4,3,2),p=reshape(repeat([:U,:V,:W,:O],outer=[6]),4,3,2))
    @fact mapslices(x->larr(:l=>[1,2,3,4],:p=>[:U,:V,:W,:O], axis1=[:m,:n,:p,:q]), darr(a=rand(3,5,2)), [2]) --> larr(l=reshape(repeat(collect(1:4),outer=[6]),4,3,2),p=reshape(repeat([:U,:V,:W,:O],outer=[6]),4,3,2),axis1=[:m,:n,:p,:q])
    @fact mapslices(x->[11,21,31,41], darr(a=rand(3,5,2)), [2]) --> reshape(repeat(nalift([11,21,31,41]),outer=[6]),4,3,2)
    @fact mapslices(x->x[:a][1,1].value==1 ? LDict(:a=>Nullable(1)) : LDict(:b=>nalift([1 2 3])), darr(a=[1 2 3;4 5 6]),[1]) --> nalift(DataCubes.simplify_array(Any[LDict(:a=>1), LDict(:b=>[1 2 3]), LDict(:b=>[1 2 3])]))
    @fact (i=0;mapslices(x->x[:a][1,1].value<=2 ? (i+=1;LDict(:a=>i)) : LDict(:b=>-1), darr(a=[1 2 3;4 5 6]),[1])) --> nalift(DataCubes.simplify_array(Any[LDict(:a=>1), LDict(:a=>2), LDict(:b=>-1)]))
    @fact (i=0;mapslices(x->x[:a][1,1].value<=2 ? (i+=1;LDict(:a=>i)) : LDict(:a=>-1), darr(a=[1 2 3;4 5 6]),[1])) --> darr(a=[1,2,-1])
    @fact mapslices(x->[1], darr(a=Int[]), [1])  --> isnull
    @fact size(mapslices(x->[1], darr(a=rand(0,5,3)), [2])) --> (0,3)
    @fact mapslices(x->msum(x,1,2), darr(a=reshape(1:24,2,3,4)),2,1) --> darr(a=reshape([1,3,6,10,15,21,7,15,24,34,45,57,13,27,42,58,75,93,19,39,60,82,105,129],2,3,4))
    @fact mapslices(x->msum(x,2,1), 3.0*darr(a=reshape(1:24,2,3,4)),[2,1]) --> 3.0*darr(a=reshape([1,3,6,10,15,21,7,15,24,34,45,57,13,27,42,58,75,93,19,39,60,82,105,129],2,3,4))

    #@fact sub(darr(a=1.0*[1 2 3;4 5 6]),2:2,1:2) --> darr(a=1.0*[4 5])
    #@fact sub(darr(a=1.0*[1 2 3;4 5 6]),(2:2,1:2)) --> darr(a=1.0*[4 5])
    @fact view(darr(a=1.0*[1 2 3;4 5 6]),2:2,1:2) --> darr(a=1.0*[4 5])
    @fact view(darr(a=1.0*[1 2 3;4 5 6]),(2:2,1:2)) --> darr(a=1.0*[4 5])
    @fact typeof(repeat(darr(a=1.0*[1 2 3;4 5 6]), inner=[2,1], outer=[1,3])) --> DataCubes.DictArray{Symbol,2,DataCubes.AbstractArrayWrapper{Nullable{Float64},2,DataCubes.FloatNAArray{Float64,2,Array{Float64,2}}},Nullable{Float64}}
    @fact reshape(darr(a=1.0*[1 3 5;2 4 6]), 1, 6) --> darr(a=1.0*[1 2 3 4 5 6])
    @fact typeof(reshape(darr(a=1.0*[1 3 5;2 4 6]), 1, 6)) --> DictArray{Symbol,2,DataCubes.AbstractArrayWrapper{Nullable{Float64},2,DataCubes.FloatNAArray{Float64,2,Array{Float64,2}}},Nullable{Float64}}
    @fact (1.0*nalift([1 2 3]))[CartesianIndex((1,2))].value --> 2.0
    @fact (a=1.0*nalift([1 2 3]);a[1,2]=5.0;a) --> 1.0*nalift([1 5 3])
    @fact (a=1.0*nalift([1 2 3]);a[2]=5.0;a) --> 1.0*nalift([1 5 3])
    @fact (a=1.0*nalift([1 2 3]);a[1:2]=[5.0 1.0];a) --> 1.0*nalift([5 1 3])
    @fact repeat(darr(a=1.0*[1 2;3 4]), inner=[2,1], outer=[1,3]) --> darr(a=1.0*[1 2 1 2 1 2;1 2 1 2 1 2;3 4 3 4 3 4;3 4 3 4 3 4])
    @fact repeat(darr(a=enumeration([1 2;3 4])), inner=[2,1], outer=[1,3]) --> darr(a=enumeration([1 2 1 2 1 2;1 2 1 2 1 2;3 4 3 4 3 4;3 4 3 4 3 4]))
    @fact repeat(darr(a=1.0*nalift([1 2;3 4])), inner=[2,1], outer=[1,3]) --> darr(a=1.0*[1 2 1 2 1 2;1 2 1 2 1 2;3 4 3 4 3 4;3 4 3 4 3 4])
    @fact repmat(darr(a=enumeration([1 2 3;4 5 6])), 2) --> darr(a=[1 2 3;4 5 6;1 2 3;4 5 6])
    @fact repmat(darr(a=1.0*nalift([1 2 3;4 5 6])), 2) --> darr(a=1.0*[1 2 3;4 5 6;1 2 3;4 5 6])
    @fact repmat(darr(a=1.0*nalift([1 2 3;4 5 6])), 2) --> darr(a=1.0*[1 2 3;4 5 6;1 2 3;4 5 6])
    @fact typeof(repmat(darr(a=enumeration([1 2;3 4])), 2, 2).data[:a]) --> DataCubes.AbstractArrayWrapper{Nullable{Int64},2,DataCubes.EnumerationArray{Int64,2,DataCubes.AbstractArrayWrapper{Int64,2,Array{Int64,2}},Int64}}
    @fact repmat(darr(a=enumeration([1 2;3 4])), 2, 2) --> darr(a=[1 2 1 2;3 4 3 4;1 2 1 2;3 4 3 4])
    @fact repmat(darr(a=1.0*nalift([1 2;3 4])), 2, 2) --> darr(a=1.0*[1 2 1 2;3 4 3 4;1 2 1 2;3 4 3 4])
    @fact map(x->LDict(:c=>x[:a]), @darr(a=[1,2,3],b=[4,5,6])) --> @darr(c=[1,2,3])
    @fact map(x->x[:a], @darr(a=[1,2,3],b=[4,5,6])) --> nalift([1,2,3])
    @fact reverse(@darr(a=[1 2 3;4 5 6])) --> @darr(a=[4 5 6;1 2 3])
    @fact reverse(@darr(a=[1 2 3;4 5 6]),[1]) --> @darr(a=[4 5 6;1 2 3])
    @fact reverse(@darr(a=[1 2 3;4 5 6]),1:2) --> @darr(a=[6 5 4;3 2 1])
    @fact reverse(@darr(a=[1 2 3;4 5 6]),[2,1]) --> @darr(a=[6 5 4;3 2 1])
    @fact reverse(@darr(a=[1 2 3;4 5 6]),2:-1:1) --> @darr(a=[6 5 4;3 2 1])
    @fact flipdim(@darr(a=[1 2 3;4 5 6]), 1, 2) --> @darr(a=[6 5 4;3 2 1])
    @fact flipdim(@darr(a=[1 2 3;4 5 6]), 1) --> @darr(a=[4 5 6;1 2 3])
    @fact larr(darr(a=1:10), axis1=11:20) --> larr(larr(a=1:10), axis1=11:20)
    @fact larr(darr(a=1:10), axis1=11:20) --> larr(a=1:10, axis1=11:20)
    @fact fill(LDict('a'=>1, 'b'=>10), 10, 20) --> darr('a'=>fill(1,10,20), 'b'=>fill(10,10,20))
    @fact fill(LDict('a'=>1, 'b'=>10), (10, 20)) --> darr('a'=>fill(1,10,20), 'b'=>fill(10,10,20))
    @fact reshape(darr(a=[1 2 3;4 5 6],b=[10 11 12;13 14 15]), 3, 2) --> darr(a=[1 5;4 3;2 6],b=[10 14;13 12;11 15])
    @fact reshape(darr(a=enumeration([1 2 3;4 5 6]),b=enumeration([10 11 12;13 14 15])), 3, 2) --> darr(a=[1 5;4 3;2 6],b=[10 14;13 12;11 15])
    @fact reducedim((x,y)->LDict(:a=>Nullable(x[:a].value+y[:a].value)), darr(a=[1 2 3;4 5 6]),[1],LDict(:a=>Nullable(0))) --> darr(a=[5,7,9])
    @fact reducedim((x,y)->x+y[:a].value, darr(a=[1 2 3;4 5 6]),[1],0) --> nalift([5,7,9])
    @fact reducedim((x,y)->x+y[:a].value, darr(a=[1 2 3;4 5 6]),[2],0) --> nalift([6,15])
    @fact reducedim((x,y)->x+y[:a].value, darr(a=[1 2 3;4 5 6]),[1,2],0).value --> 21
    @fact reorder(darr('x'=>1:10, 3=>11:20), 3) --> darr(3=>11:20, 'x'=>1:10)
    @fact reorder(darr(c1=1:10,c2=11:20),:c2,:c1) --> reorder(darr(c1=1:10,c2=11:20),:c2)
    @fact reorder(darr(c1=1:10,c2=11:20),:c2,:c1) --> darr(c2=11:20,c1=1:10)
    @fact rename(darr('x'=>1:10, 3=>11:20), 1) --> darr(1=>1:10, 3=>11:20)
    @fact rename(darr(c1=1:10,c2=11:20),:a,:b) --> darr(a=1:10,b=11:20)
    @fact cat(1, darr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), darr(b=['x' 'y' 'z'], d=[:m :n :p])) --> reshape(@darr(a=[1,4,NA,2,5,NA,3,6,NA],b=['a','d','x','b','e','y','c','f','z'],d=[NA,NA,:m,NA,NA,:n,NA,NA,:p]), 3, 3)
    @fact cat(2, darr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), darr(b=['x','y'], d=[:m,:n])) --> reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f','x','y'], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    @fact cat(2, darr(a=[1 2 3;4 5 6], b=['a' 'b' 'c';'d' 'e' 'f']), darr(b=[10,11], d=[:m,:n])) --> reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    # TODO: decide whether to allow type promotion from Nullabe{T} to Nullable{U} if T can be promoted to U.
    # currently, do not do that.
    # @fact cat(1,darr(k=[1 2 3]),darr(k=[4.0 5.0 6.0])) --> darr(k=[1.0 2.0 3.0;4.0 5.0 6.0])
    @fact cat(1,darr(k=[1.0 2.0 3.0]),darr(k=[4.0 5.0 6.0])) --> darr(k=[1.0 2.0 3.0;4.0 5.0 6.0])
    @fact vcat(darr(k=[1.0 2.0 3.0]),darr(k=[4.0 5.0 6.0])) --> darr(k=[1.0 2.0 3.0;4.0 5.0 6.0])
    @fact hcat(darr(k=[1.0 2.0 3.0]),darr(k=[4.0 5.0 6.0])) --> darr(k=[1.0 2.0 3.0 4.0 5.0 6.0])
    @fact cat(1, darr(a=enumeration([1 2 3;4 5 6]), b=enumeration(['a' 'b' 'c';'d' 'e' 'f'],['x','y','z','a','b','c','d','e','f'])), darr(b=enumeration(['x' 'y' 'z'],['x','y','z','a','b','c','d','e','f']), d=[:m :n :p])) --> reshape(@darr(a=[1,4,NA,2,5,NA,3,6,NA],b=['a','d','x','b','e','y','c','f','z'],d=[NA,NA,:m,NA,NA,:n,NA,NA,:p]), 3, 3)
    @fact vcat(darr(a=enumeration([1 2 3;4 5 6]), b=enumeration(['a' 'b' 'c';'d' 'e' 'f'],['x','y','z','a','b','c','d','e','f'])), darr(b=enumeration(['x' 'y' 'z'],['x','y','z','a','b','c','d','e','f']), d=[:m :n :p])) --> reshape(@darr(a=[1,4,NA,2,5,NA,3,6,NA],b=['a','d','x','b','e','y','c','f','z'],d=[NA,NA,:m,NA,NA,:n,NA,NA,:p]), 3, 3)
    @fact cat(2, darr(a=[1 2 3;4 5 6], b=enumeration(['a' 'b' 'c';'d' 'e' 'f'],['a','b','c','d','e','f'])), darr(b=enumeration(['a','b'],['a','b','c','d','e','f']), d=[:m,:n])) --> reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f','a','b'], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    @fact cat(2, darr(a=enumeration([1 2 3;4 5 6]), b=['a' 'b' 'c';'d' 'e' 'f']), darr(b=[10,11], d=[:m,:n])) --> reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    @fact hcat(darr(a=enumeration([1 2 3;4 5 6]), b=['a' 'b' 'c';'d' 'e' 'f']), darr(b=[10,11], d=[:m,:n])) --> reshape(@darr(a=[1,4,2,5,3,6,NA,NA], b=['a','d','b','e','c','f',10,11], d=[NA,NA,NA,NA,NA,NA,:m,:n]), 2, 4)
    @fact reshape(darr(a=enumeration(11:16)),2,3) --> darr(a=[11 13 15;12 14 16])
    @fact reshape(darr(a=enumeration(11:16)),(2,3)) --> darr(a=[11 13 15;12 14 16])
    @fact mapna((x,y)->(x,y),nalift([1 2 3;4 5 6]),nalift(1.0*[2 3 4;5 6 7])) --> nalift([(1,2.0) (2,3.0) (3,4.0);(4,5.0) (5,6.0) (6,7.0)])
    @fact show(darr(a=rand(2))) --> nothing
    @fact show(darr(a=rand(2,3))) --> nothing
    @fact show(darr(a=rand(2,3,2))) --> nothing
    @fact show(darr(a=rand(2,3,2,2))) --> nothing

    context("show tests") do
      @fact show(darr(a=view([1,2],1))) --> nothing
      @fact show(darr(a=[])) --> nothing
      @fact show(darr(a=rand(2))) --> nothing
      @fact show(darr(a=rand(2,3))) --> nothing
      @fact show(darr(a=rand(2,3,2))) --> nothing
      @fact show(darr(a=rand(2,3,2,2))) --> nothing
      @fact show(darr(a=rand(2,3,2,2))) --> nothing
      @fact (dcube.set_showsize!(5,5);show(darr(a=rand(10,10)))) --> nothing
      @fact (dcube.set_showheight!(3);show(darr(a=rand(10,10)))) --> nothing
      @fact (dcube.set_showwidth!(3);show(darr(a=rand(10,10)))) --> nothing
      @fact (dcube.set_default_showsize!();nothing) --> nothing
      @fact (dcube.set_showalongrow!(false);show(darr(a=rand(3),b=rand(3),c=fill(:X,3)))) --> nothing
      @fact (dcube.set_showalongrow!(true);show(darr(a=rand(3),b=rand(3),c=fill(:X,3)))) --> nothing
      @fact (dcube.set_showalongrow!(false);show(darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_showalongrow!(true);show(darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_format_string!(Float64, "%0.2f");show(darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_format_string!(Float64, "%0.8g");show(darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_dispsize!(5,5);show(STDOUT,MIME("text/html"),darr(a=rand(10,10)))) --> nothing
      @fact (dcube.set_dispheight!(3);show(STDOUT,MIME("text/html"),darr(a=rand(10,10)))) --> nothing
      @fact (dcube.set_dispwidth!(3);show(STDOUT,MIME("text/html"),darr(a=rand(10,10)))) --> nothing
      @fact (dcube.set_dispwidth!(3);show(STDOUT,MIME("text/html"),darr(a=rand(2,3,4)))) --> nothing
      @fact (dcube.set_dispwidth!(3);show(STDOUT,MIME("text/html"),darr(a=view([1,2],1)))) --> nothing
      @fact (dcube.set_default_dispsize!();nothing) --> nothing
      @fact (dcube.set_dispalongrow!(false);show(STDOUT,MIME("text/html"),darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_dispalongrow!(true);show(STDOUT,MIME("text/html"),darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_format_string!(Float64, "%0.2f");show(STDOUT,MIME("text/html"),darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
      @fact (dcube.set_format_string!(Float64, "%0.8g");show(STDOUT,MIME("text/html"),darr(a=rand(3,5),b=rand(3,5),c=fill(:X,3,5)))) --> nothing
    end
  end
end

end