text:
Lar = CuboidGrids
using DataStructures
using Distributed

"""
	INSR(f::Function)(seq::Array{Any,1})::Any

FL primitive combinator to transform a binary function to an n-ary one.
```
julia> mod1D = Lar.grid(repeat([.1,-.1],outer=5)...)
([0.0 0.1 … 0.9 1.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [9, 10]])

julia> using ViewerGL; GL = ViewerGL

julia> GL.VIEW([ GL.GLFrame2, GL.GLGrid(mod1D..., GL.COLORS[1],1) ])

julia> mod3D = Lar.INSR(Lar.larModelProduct)([mod1D,mod1D,mod1D])
([0.0 0.0 … 1.0 1.0; 0.0 0.0 … 1.0 1.0; 0.0 0.1 … 0.9 1.0],
Array{Int64,1}[[1, 2, 12, 13, 122, 123, 133, 134], [3, 4, 14, 15, 124, 125, 135, 136],
… [1063, 1064, 1074, 1075, 1184, 1185, 1195, 1196], [1065, 1066, 1076, 1077, 1186, 1187, 1197, 1198]])

julia> GL.VIEW([ GL.GLFrame2, GL.GLPol(mod3D..., GL.COLORS[1],1) ])
```
"""
function INSR(f)
	function INSR0(seq)
		len = length(seq)
		res = seq[end]
		for i in range(len-2,step=-1,stop=0)
			res = f([seq[i+1], res])
		end
		return res
	end
	return INSR0
end

"""
	grid(sequence::Array{Number,1})::Lar.LAR

Generate a 1D LAR model.

*Geometry* is stored as 1D `Points`, and *Topology* is stored as 1D `Cells`.
`q()` and `q()()` are used as alias function name.
```julia
julia> model1D = Lar.grid(1,-1,1,-1,1,-1,1,-1,1,-1)
([0.0 1.0 … 9.0 10.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [9, 10]])

julia> model1D[1]
1×11 Array{Float64,2}:
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  9.0  10.0

julia> model1D[2]
5-element Array{Array{Int64,1},1}:
 [1, 2]
 [3, 4]
 [5, 6]
 [7, 8]
 [9, 10]

 julia> mod = Lar.grid(repeat([.1,-.1],outer=5)...)
 ([0.0 0.1 … 0.9 1.0], Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [9, 10]])
```
"""
function grid(sequence...)
	sequence = collect(sequence)
	cursor,points,hulls= (0,[[0.]],[])
	for value in sequence
		points = append!(points, [[cursor + abs(value)]])
		if value>=0
			append!(hulls,[[length(points)-1,length(points)]])
		end
	  cursor += abs(value)
	end
	V = convert(Lar.Points, [p[1] for p in points]')
	EV = convert(Lar.Cells,hulls)
	return V,EV
end
const q = grid

"""
	qn(n::Int)(sequence::Array{T,1})::Lar.LAR  where T <: Real

Alias of `grid` function, with repetition parameter `n`.
```
julia> Lar.qn(3)([1.5,-2,0.5])
([0.0 1.5 … 11.5 12.0], Array{Int64,1}[[1, 2], [3, 4], [4, 5], [6, 7], [7, 8], [9, 10]])
```
"""
function qn(n::Int)
	function qn0(sequence::Array{T,1})::Lar.LAR  where T <: Real
		sequence = collect(sequence)
		return Lar.grid(repeat(sequence,outer=n)...)
	end
	return qn0
end

"""
	grid_0(n::Int)::Array{Int64,2}

Generate a *uniform 0D cellular complex*.
The `grid_0` function generates a 0-dimensional uniform complex embedding ``n+1`` equally-spaced  0-cells (at *unit interval* boundaries). It returns by columns the cells of this 0-complex as `Array{Int64,2}.

#	Example
```julia
julia> grid_0(10)
# output
1×11 Array{Int64,2}:
 0  1  2  3  4  5  6  7  8  9  10
```
"""
function grid_0(n::Int)::Array{Int64,2}
    return hcat([[i] for i in range(0, length=n+1)]...)
end

"""
	grid_1(n::Int)::Array{Int64,2}

Generate a *uniform 1D cellular complex*.
The `grid_1` function generates a 0-dimensional uniform complex embedding ``n+1`` equally-spaced  1-cells (*unit intervals*). It returns by columns the cells of this 1-complex as `Array{Int64,2}`.

#	Example
```julia
julia> grid_1(10)
# output
2×10 Array{Int64,2}:
 0  1  2  3  4  5  6  7  8   9
 1  2  3  4  5  6  7  8  9  10
```
"""
function grid_1(n)
    return hcat([[i,i+1] for i in range(0, length=n)]...)
end

"""
	larGrid(n::Int)(d::Int)::Array{Int64,2}

Generate either a *uniform 0D cellular complex* or a *uniform 1D cellular complex*.
A `larGrid` function is given to generate the LAR representation of the cells of either a 0- or a 1-dimensional complex, depending on the value of the `d` parameter, to take values in the set ``{0,1}``, and providing the *order* of the output complex.

#	Example
```julia
julia> larGrid(10)(0)
# output
1×11 Array{Int64,2}:
 0  1  2  3  4  5  6  7  8  9  10

julia> larGrid(10)(1)
# output
2×10 Array{Int64,2}:
 0  1  2  3  4  5  6  7  8   9
 1  2  3  4  5  6  7  8  9  10
```
"""
function larGrid(n::Int)
    function larGrid1(d::Int)::Array{Int64,2}
        if d==0
         return grid_0(n)
        elseif d==1
         return grid_1(n)
        end
    end
    return larGrid1
end

"""
	cart(args::Array{Array{Any,1},1})::Array{Tuple,1}

Cartesian product of collections given in the
 unary `Array` argument. Return an `Array` of `Tuple`. The number ``n`` of output `Tuple` is equal to the *product of sizes* of input `args`

#	Example
```julia
julia> cart([[1,2,3],["a","b"],[11,12]])
# output
12-element Array{Tuple{Any,Any,Any},1}:
 (1, "a", 11)
 (1, "a", 12)
 (1, "b", 11)
 (1, "b", 12)
 (2, "a", 11)
 (2, "a", 12)
 (2, "b", 11)
 (2, "b", 12)
 (3, "a", 11)
 (3, "a", 12)
 (3, "b", 11)
 (3, "b", 12)
```
"""
function cart(args)::Array{Tuple,1}
   return sort(vcat(collect(Iterators.product(args...))...))
end

"""
	larVertProd(vertLists::Array{Points,1})::Points

Generate the integer *coordinates of vertices* (0-cells) of a *multidimensional grid*.
*Grid n-vertices* are produced by the `larVertProd` function, via Cartesian product of vertices of ``n`` 0-dimensional arguments (vertex arrays in `vertLists`), orderly corresponding to ``x_1, x_2, ..., x_n`` coordinates in the output points ``(x_1, x_2,...,x_n)`` in ``R^n``.

#	Example

```julia
julia> larVertProd([ larGrid(3)(0), larGrid(4)(0) ])
# output
2×20 Array{Int64,2}:
 0  0  0  0  0  1  1  1  1  1  2  2  2  2  2  3  3  3  3  3
 0  1  2  3  4  0  1  2  3  4  0  1  2  3  4  0  1  2  3  4
```
"""
function larVertProd(vertLists::Array{Array{Int64,2},1})::Array{Int64,2}
   coords = [[x[1] for x in v] for v in Lar.cart(vertLists)]
   return sortslices(hcat(coords...), dims=2)
end
function larVertProd(vertLists::Array{Array{Float64,2},1})::Array{Float64,2}
   coords = [[x[1] for x in v] for v in Lar.cart(vertLists)]
   return sortslices(hcat(coords...), dims=2)
end

"""
	index2addr(shape::Array{Int64,1})(multiIndex)::Int

*Multi-index to address* transformation. Multi-index is a *generalization* of the concept of an integer index to an *ordered tuple of indices*.
The second-order utility function `index2addr`  transforms a `shape` list for a *multidimensional array* into a function that, when applied to a *multindex array*, i.e. to a list of integer `Tuple` within the `shape`'s bounds, returns the *integer addresses* of the corresponding array components within the *linear storage* of the multidimensional array.

#	Example
Notice that in the example below, there are ``3 x 6`` different *multi-index values* for the variable `index`, generated by `cart([ 0:2, 0:5 ])`.

```julia
julia> [index2addr([3,6])(collect(index)) for index in cart([ 0:2, 0:5 ])]'
# output
1×18 RowVector{Int64,Array{Int64,1}}:
 1  2  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18

julia> index2addr([3,6])([0,0])
# output
1

julia> index2addr([3,6])([2,5])
# output
18
```
"""
function index2addr( shape::Array{Int64,2} )
    n = length(shape)
    theShape = append!(shape[2:end],1)
    weights = [prod(theShape[k:end]) for k in range(1, length=n)]

    function index2addr0( multiIndex::Array{Int,1} )::Int
        return dot(collect(multiIndex), weights) + 1
    end

    return index2addr0
end

"""
	index2addr(shape::Array{Int64,2})(multiIndex::Array{Int,1})::Int

Multi-index to address transformation. Partial function allowing for using both horizontal and vertical vectors for `shape` parameter. Notice that multi-indices are used here as *coordinates of grid points*, hence they start from tuples of zeros. Accordingly, the translation formula for multi-index to address transformation is *0-based*.
"""
function index2addr(shape::Array{Int64,1})
   index2addr(hcat(shape...))
end

"""
	larCellProd(cellLists::Array{Cells,1})::Cells

Generation of *grid cells* by *Cartesian product* of 0/1-complexes.
The *output complex* is generated by the product of *any number* of either 0- or 1-dimensional cell complexes. The product of ``d`` 1-complexes generates *solid ``d``-cells*, while the product of ``n`` 0-complexes and ``d-n`` 1-complexes (``n < d``) generates *non-solid ``(d-n)``-cells*, properly embedded in ``d``-space, i.e. with vertices having ``d`` coordinates.

# Examples

To understand the *generation of cuboidal grids* from products of 0- or 1-dimensional complexes, below we show a simple example of 2D grids embedded in ``R^3``.
In particular, `v1 = [0. 1. 2. 3.]` and `v0 = [0. 1. 2.]` are two 2-arrays of 1D vertices, `c1 = [[0,1],[1,2],[2,3]]` and `c0 = [[0],[1],[2]]` are the LAR representation of one *``1``-complex* and one *``0``-complex*, respectively. The solid 2-complex named `grid2D` is generated in 2D as follows:

```julia
julia> v1 = [0. 1. 2. 3.]
1×4 Array{Float64,2}:
 0.0  1.0  2.0  3.0

julia> c1 = [[0,1],[1,2],[2,3]]
3-element Array{Array{Int64,1},1}:
 [0, 1]
 [1, 2]
 [2, 3]

julia> grid2D = larVertProd([v1,v1]),larCellProd([c1,c1])
([0.0 0.0 … 3.0 3.0; 0.0 1.0 … 2.0 3.0], Array{Int64,1}[[1, 2, 5, 6], [2, 3, 6, 7], [3, 4, 7, 8], [5, 6, 9, 10], [6, 7, 10, 11], [7, 8, 11, 12], [9, 10, 13, 14], [10, 11, 14, 15], [11, 12, 15, 16]])
```
whereas a *non-solid* ``2``-complex in ``3D`` is generated as:

```julia
julia> v1, c1 = [0. 1. 2. 3.],[[0,1],[1,2],[2,3]]
([0.0 1.0 2.0 3.0], Array{Int64,1}[[0, 1], [1, 2], [2, 3]])

julia> v0, c0 = [0. 1. 2.], [[0],[1],[2]]
([0.0 1.0 2.0], Array{Int64,1}[[0], [1], [2]])

julia> vertGrid = larVertProd([v1, v1, v0])
3×48 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  …  3.0  3.0  3.0  3.0  3.0  3.0  3.0  3.0  3.0
 0.0  0.0  0.0  1.0  1.0  1.0  …  1.0  1.0  1.0  2.0  2.0  2.0  3.0  3.0  3.0
 0.0  1.0  2.0  0.0  1.0  2.0  …  0.0  1.0  2.0  0.0  1.0  2.0  0.0  1.0  2.0

julia> cellGrid = larCellProd([c1, c1, c0])
27-element Array{Array{Int64,1},1}:
 [1, 4, 13, 16]
 [2, 5, 14, 17]
 ...  ... ...
 [32, 35, 44, 47]
 [33, 36, 45, 48]

julia> grid3D = vertGrid,cellGrid
([0.0 0.0 … 3.0 3.0; 0.0 0.0 … 3.0 3.0; 0.0 1.0 … 1.0 2.0], Array{Int64,1}[[1, 4, 13, 16], [2, 5, 14, 17], … [32, 35, 44, 47], [33, 36, 45, 48]])

julia> using Plasm

julia> Plasm.view(grid3D)
```
"""
function larCellProd(cellLists::Array{Cells,1})::Cells
   shapes = [length(item) for item in cellLists]
   subscripts = cart([collect(range(0, length=shape)) for shape in shapes])
   indices = [collect(tuple) .+ 1 for tuple in subscripts]

   jointCells = [cart([cells[k] for (k,cells) in zip(index,cellLists)])
   				for index in indices]
   convertIt = index2addr([ (length(cellLists[k][1]) > 1) ? shape .+ 1 : shape
      for (k,shape) in enumerate(shapes) ])
   [vcat(map(convertIt, map(collect,jointCells[j]))...) for j in 1:length(jointCells)]
   #[vcat(pmap(convertIt, pmap(collect,jointCells[j]))...) for j in 1:length(jointCells)] la parlallelizzazione appesantisce l'esecuzione(?)

end

"""
	filterByOrder( n::Int )Array{Array{Array{Int8,1},1},1}

Filter the `array` of codes  (`Boolean` `String`) of *``n`` bits* depending on their integer value (*order*).

# Example

```julia
julia> filterByOrder(3)
# output
4-element Array{Array{Array{Int8,1},1},1}:
 Array{Int8,1}[Int8[0, 0, 0]]
 Array{Int8,1}[Int8[0, 0, 1], Int8[0, 1, 0], Int8[1, 0, 0]]
 Array{Int8,1}[Int8[0, 1, 1], Int8[1, 0, 1], Int8[1, 1, 0]]
 Array{Int8,1}[Int8[1, 1, 1]]
```"""
function filterByOrder(n::Int)Array{Array{Array{Int8,1},1},1}
   terms = [[parse(Int8,bit) for bit in collect(term)] for term in Lar.binaryRange(n)]
   return [[term for term in terms if sum(term) == k] for k in 0:n]
end

"""
	larGridSkeleton( shape::Array{Int,1} )( d::Int )::Cells

Produce the `d`-dimensional skeleton (set of `d`-cells) of a cuboidal grid of given `shape`.

# Example

A `shape=[1,1,1]` parameter refers to a *grid* with a single step on the three axes, i.e. to a single *3D unit cube*. Below all *skeletons* of such simplest grid are generated.

```julia
julia> Lar.larGridSkeleton([1,1,1])(0)
# output
8-element Array{Array{Int64,1},1}:
[[1], [2], [3], [4], [5], [6], [7], [8]]

julia> Lar.larGridSkeleton([1,1,1])(1)
# output
12-element Array{Array{Int64,1},1}:
[[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> Lar.larGridSkeleton([1,1,1])(2)
# output
6-element Array{Array{Int64,1},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]]

julia> Lar.larGridSkeleton([1,1,1])(3)
# output
1-element Array{Array{Int64,1},1}:
 [1, 2, 3, 4, 5, 6, 7, 8]
```
"""
function larGridSkeleton(shape)
    n = length(shape)
    function larGridSkeleton0( d::Int )::Cells

    	@assert d<=n

        components = filterByOrder(n)[d .+ 1]
        apply(fun,a) = fun(a)
		componentCellLists = [ [pmap(f,x) for (f,x) in zip( [larGrid(dim)
        for dim in shape], convert(Array{Int64,1},component) ) ]
				for component in components ]
        colList(arr) = [arr[:,k]  for k in 1:size(arr,2)]
        out = [ larCellProd(pmap(colList,cellLists)) for cellLists in componentCellLists ]
        return vcat(out...)
    end
    return larGridSkeleton0
end

function larGridSkeletonT(shape)
    n = length(shape)
    l = Threads.Lock()
    Threads.@spawn function larGridSkeleton0T( d::Int )::Cells
    	@assert d<=n
        components = filterByOrder(n)[d .+ 1]
		componentCellLists = [ [pmap(f,x)  for (f,x) in  zip( [larGrid(dim)
        for dim in shape], convert(Array{Int64,1},component) ) ]
				for component in components ]
        colList(arr) = [arr[:,k]  for k in 1:size(arr,2)]
        out = [ larCellProd(pmap(colList,cellLists)) for cellLists in componentCellLists ]
        return vcat(out...)
    end
    return larGridSkeleton0T
end

"""
	larImageVerts(shape::Array{Int,1})::Array{Int64,2}

Linearize the *grid of integer vertices*, given the `shape` of a *cuboidal grid* (typically an *image*).

# Examples

```julia
julia> larImageVerts([1024,1024])
# output
2×1050625 Array{Int64,2}:
 0  0  0  0  0  0  0  0  0  0   0   0 … 1024  1024  1024  1024  1024  1024  1024  1024
 0  1  2  3  4  5  6  7  8  9  10  11 … 1017  1018  1019  1020  1021  1022  1023  1024

julia> larImageVerts([1,1,1])
# output
3×8 Array{Int64,2}:
 0  0  0  0  1  1  1  1
 0  0  1  1  0  0  1  1
 0  1  0  1  0  1  0  1
```
"""
function larImageVerts( shape::Array{Int,1} )::Array{Int64,2}
   vertexDomain(n) = hcat([k for k in 0:n-1]...)
   vertLists = [vertexDomain(k+1) for k in shape]
   #vertList = [pmap(vertexDomain, shape)+1] qui non si puo mettere pmap
   vertGrid = larVertProd(vertLists)
   return vertGrid
end

#°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°°

function vertexDomainAux(n::Int)::Array{Array{Float64, 2}, 1}
    ret = Array
    for k in n-1
        ret = hcat(k...)
        end
    return ret
    end

function larImageVertsT( shape::Array{Int,1} )::Array{Int64,2}
    vertLists = [vertexDomainAux(k+1) for k in shape]
    vertGrid = larVertProd(vertLists)
    return vertGrid
 end

"""
	cuboidGrid( shape, filled=false )::Union( Cells, Array{Cells,1} )

Multi-dimensional generator function.
Generate either a solid *``d``-grid* of unit *``d``-cuboids* in ``d``-dimensional space, or the array of ``p``-skeletons (``0 <=p<= d``), depending on the Boolean variable `filled`. ``0``-cuboids are points, ``1``-cuboids are segments, , ``2``-cuboids are squares,  ``3``-cuboids are cubes, etc. The `shape=[a,b,c]` value determines the number ``a x b x c`` of ``d``-cells. Notice that `d = length(shape)`
"""
function cuboidGrid( shape, filled=false )
   vertGrid = larImageVerts(shape)
   gridMap = larGridSkeleton(shape)
   #fetch(vertGrid)
   if ! filled
      cells = gridMap(length(shape))
   else
      skeletonIds = 0:length(shape)
      #cells = [ gridMap(id) for id in skeletonIds ]
      cells = [ pmap(gridMap, skeletonIds )]
   end
   return convert(Array{Float64,2},vertGrid), cells
end

function cuboidGridT( shape, filled=false )
    vertGrid = larImageVertsT(shape)
    gridMap = larGridSkeletonT(shape)

    #sul ramo di else non ci arrivo mai oggettivo fattuale
    cells = gridMap(length(shape))

    return convert(Array{Float64,2},vertGrid), cells
 end

"""
	larModelProduct

The `larModelProduct` function takes as input a pair of *LAR models* and returns the model of their *Cartesian product*. Since LAR type is a pair ``(geometry,topology)``, the second element of output is the *topological product* of the input topologies.

# Example

Data preparation follows.

```julia
julia> geom_0,topol_0 = [0. 1. 2. 3. 4.],[[1],[2],[3],[4],[5]]
([0.0 1.0 … 3.0 4.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])

julia> geom_1,topol_1 = [0. 1. 2.], [[1,2],[2,3]]
([0.0 1.0 2.0], Array{Int64,1}[[1, 2], [2, 3]])

julia> mod_0 = (geom_0,topol_0)
([0.0 1.0 … 3.0 4.0], Array{Int64,1}[[1, 2], [2, 3], [3, 4], [4, 5]])

julia> mod_1 = (geom_1,topol_1)
([0.0 1.0 2.0], Array{Int64,1}[[1, 2], [2, 3]])
```
Generation of a 2D `squares` model, with 8 *two-dimensional cells*.

```julia
julia> squares = larModelProduct(mod_1,mod_1)
([0.0 0.0 … 4.0 4.0; 0.0 1.0 … 1.0 2.0], Array{Int64,1}[[1, 2, 4, 5], [2, 3, 5, 6], [4, 5, 7, 8], [5, 6, 8, 9], [7, 8, 10, 11], [8, 9, 11, 12], [10, 11, 13, 14], [11, 12, 14, 15]])

julia> squares[1]
2×15 Array{Float64,2}:
 0.0  0.0  0.0  1.0  1.0  1.0  2.0  2.0  2.0  3.0  3.0  3.0  4.0  4.0  4.0
 0.0  1.0  2.0  0.0  1.0  2.0  0.0  1.0  2.0  0.0  1.0  2.0  0.0  1.0  2.0

julia> squares[2]
8-element Array{Array{Int64,1},1}:
[[1,2,4,5], [2,3,5,6], [4,5,7,8], [5,6,8,9], [7,8,10,11], [8,9,11,12], [10,11,13,14], [11,12,14,15]]
```

Generation of a 3D `cubes` model, with 32 *three-dimensional cells*.

```julia
julia> cubes = larModelProduct(squares,mod_0)
([0.0 0.0 … 4.0 4.0; 0.0 0.0 … 2.0 2.0; 0.0 1.0 … 3.0 4.0], Array{Int64,1}[[1, 2, 6, 7, 16, 17, 21, 22], [2, 3, 7, 8, 17, 18, 22, 23], [3, 4, 8, 9, 18, 19, 23, 24], [4, 5, 9, 10, 19, 20, 24, 25], … [53, 54, 58, 59, 68, 69, 73, 74], [54, 55, 59, 60, 69, 70, 74, 75]])
```
"""
function larModelProduct( modelOne, modelTwo )
    (V, cells1) = modelOne
    (W, cells2) = modelTwo

    vertices = DataStructures.OrderedDict();
    k = 1
    for j in 1:size(V,2)
       v = V[:,j]
        for i in 1:size(W,2)
          w = W[:,i]
            id = [v;w]
            if haskey(vertices, id) == false
                vertices[id] = k
                k = k + 1
            end
        end
    end

    cells = []
    for c1 in cells1
        for c2 in cells2
            cell = []
            for vc in c1
                for wc in c2
                    push!(cell, vertices[[V[:,vc];W[:,wc]]] )
                end
            end
            push!(cells, cell)
        end
    end

    vertexmodel = []
    for v in keys(vertices)
        push!(vertexmodel, v)
    end
    verts = hcat(vertexmodel...)
    cells = [[v for v in cell] for cell in cells]
    return (verts, cells)
end

"""
	larModelProduct(twoModels)

Further *method* associated to `larModelProduct` *function*.
"""
function larModelProduct( twoModels )
    modelOne, modelTwo = twoModels
    larModelProduct(modelOne, modelTwo)
end


text:
using ParallelAccelerator

@acc function kernel_gesummv(n::Int64, alpha::Float64, beta::Float64,
                             A::Matrix{Float64}, B::Matrix{Float64},
                             x::Vector{Float64})
    y = Array{Float64}(n)
    for i = 1:n
        t = 0.0
        y[i] = 0.0
        for j = 1:n
            t = A[i, j] * x[j] + t
            y[i] = B[i, j] * x[j] + y[i]
        end
        y[i] = alpha * t + beta * y[i]
    end
    return y
 end

@inbounds function kernel_gesummv_o(n::Int64, alpha::Float64, beta::Float64,
                             A::Matrix{Float64}, B::Matrix{Float64},
                             x::Vector{Float64})
    tic()
    y = Array{Float64}(n)
    for i = 1:n
        t = 0.0
        y[i] = 0.0
        for j = 1:n
            t = A[i, j] * x[j] + t
            y[i] = B[i, j] * x[j] + y[i]
        end
        y[i] = alpha * t + beta * y[i]
    end
    println("Julia ", toq())
    return y
 end

 n = 15000
 alpha = 1.5
 beta = 1.2

 A = rand(n, n)
 B = rand(n, n)
 x = rand(n)

res_tiramisu = kernel_gesummv(n, alpha, beta, A, B, x)
res_julia = kernel_gesummv_o(n, alpha, beta, A, B, x)

function compare(A, B)
    for i in 1:n
        if A[i] != B[i]
            println("different $((i)) : $(A[i]) != $(B[i])")
            return false
        end
    end
    return true
end

println(compare(res_tiramisu, res_julia))

ext: jl

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic
using IntervalContractors

#Preamble
setprecision(53)
setprecision(Interval, Float64)

# Using approximate checks for validaty update later?
import Base.isapprox
isapprox(x::Interval,y::Interval) = isapprox(x.lo,y.lo,atol=1E-4) && isapprox(x.hi,y.hi,atol=1E-4)

@testset "asinh_rev_test" begin
    @test asinh_rev(∅, -∞..∞)[2] == ∅
    @test isapprox(asinh_rev(Interval(0.0, Inf), -∞..∞)[2], Interval(0.0, ∞))
    @test isapprox(asinh_rev(Interval(0.0, 1.0), -∞..∞)[2], Interval(0, 1.17521))
    @test isapprox(asinh_rev(Interval(-0.5, 1.0), -∞..∞)[2], Interval(-0.521096, 1.17521))
    @test isapprox(asinh_rev(Interval(-1000.0, 1.0), -∞..∞)[2], Interval(-∞, 1.17521))
end

@testset "acosh_rev_test" begin
    @test acosh_rev(∅, -∞..∞)[2] == ∅
    @test acosh_rev(Interval(0.0, Inf), -∞..∞)[2] == Interval(1.0, ∞)
    @test isapprox(acosh_rev(Interval(0.0, 1.0), -∞..∞)[2], Interval(1.0, 1.54309))
    @test isapprox(acosh_rev(Interval(-0.5, 1.0), -∞..∞)[2], Interval(1.0, 1.54309))
    @test acosh_rev(Interval(-1000.0, 1.0), -∞..∞)[2] == Interval(1.0, ∞)
end

@testset "atanh_rev_test" begin
    @test atanh_rev(∅, -∞..∞)[2] == ∅
    @test isapprox(atanh_rev(Interval(-10.0, -1.0), -∞..∞)[2], Interval(-1, -0.761594))
    @test atanh_rev(Interval(0.0, Inf), -∞..∞)[2] == Interval(0, 1)
    @test isapprox(atanh_rev(Interval(0.0, 1.0), -∞..∞)[2], Interval(0, 0.761595))
    @test isapprox(atanh_rev(Interval(-0.5, 1.0), -∞..∞)[2], Interval(-0.462118, 0.761595))
    @test isapprox(atanh_rev(Interval(-1000.0, 1.0), -∞..∞)[2], Interval(-1, 0.761595))
    @test atanh_rev(Interval(0.0, 25.0), -∞..∞)[2] == Interval(0, 1)
    @test isapprox(atanh_rev(Interval(-1.0, 25.0), -∞..∞)[2], Interval(-0.761595, 1))
end


text:
struct NewtonRaphson{CS, AD, DT, L} <: AbstractNewtonAlgorithm{CS,AD}
    diff_type::DT
    linsolve::L
end

function NewtonRaphson(;autodiff=true,chunk_size=12,diff_type=Val{:forward},linsolve=DEFAULT_LINSOLVE)
    NewtonRaphson{chunk_size, autodiff, typeof(diff_type), typeof(linsolve)}(diff_type, linsolve)
end

mutable struct NewtonRaphsonCache{ufType, L, jType, uType, JC}
    uf::ufType
    linsolve::L
    J::jType
    du1::uType
    jac_config::JC
end

function alg_cache(alg::NewtonRaphson, f, u, p, ::Val{true})
    uf = JacobianWrapper(f,p)
    linsolve = alg.linsolve(Val{:init}, f, u)
    J = false .* u .* u'
    du1 = zero(u)
    tmp = zero(u)
    if alg_autodiff(alg)
        jac_config = ForwardDiff.JacobianConfig(uf, du1, u)
    else
        if alg.diff_type != Val{:complex}
            du2 = zero(u)
            jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg.diff_type)
        else
            jac_config = FiniteDiff.JacobianCache(Complex{eltype(tmp)}.(tmp),Complex{eltype(du1)}.(du1),nothing,alg.diff_type,eltype(u))
        end
    end
    NewtonRaphsonCache(uf, linsolve, J, du1, jac_config)
end

function alg_cache(alg::NewtonRaphson, f, u, p, ::Val{false})
    nothing
end

function perform_step(solver::NewtonImmutableSolver, alg::NewtonRaphson, ::Val{true})
    @unpack u, fu, f, p, cache = solver
    @unpack J, linsolve, du1 = cache
    calc_J!(J, solver, cache)
    # u = u - J \ fu
    linsolve(du1, J, fu, true)
    @. u = u - du1
    f(fu, u, p)
    if solver.internalnorm(solver.fu) < solver.tol
        @set! solver.force_stop = true
    end
    return solver
end

function perform_step(solver::NewtonImmutableSolver, alg::NewtonRaphson, ::Val{false})
    @unpack u, fu, f, p = solver
    J = calc_J(solver, ImmutableJacobianWrapper(f, p))
    @set! solver.u = u - J \ fu
    fu = f(solver.u, p)
    @set! solver.fu = fu
    if iszero(solver.fu) || solver.internalnorm(solver.fu) < solver.tol
        @set! solver.force_stop = true
    end
    return solver
end


text:
using DimensionalData, Test, Unitful
using DimensionalData: _slicespan, _slicemode, isrev, _bounds

@testset "order" begin
    @test indexorder(Ordered()) == ForwardIndex()
    @test arrayorder(Ordered()) == ForwardArray()
    @test relation(Ordered()) == ForwardRelation()
    @test indexorder(Unordered()) == UnorderedIndex()
    @test arrayorder(Unordered()) == ForwardArray()
    @test relation(Unordered()) == ForwardRelation()
    @test indexorder(Sampled(; order=Ordered())) == ForwardIndex()
    @test arrayorder(Sampled(; order=Ordered())) == ForwardArray()
    @test relation(Sampled(; order=Ordered())) == ForwardRelation()
end

@testset "isrev" begin
    @test isrev(ForwardIndex()) == false
    @test isrev(ForwardArray()) == false
    @test isrev(ForwardRelation()) == false
    @test isrev(ReverseIndex()) == true
    @test isrev(ReverseArray()) == true
    @test isrev(ReverseRelation()) == true
end

@testset "reverse" begin
    @test reverse(ReverseIndex()) == ForwardIndex()
    @test reverse(ForwardIndex()) == ReverseIndex()
    @test reverse(ReverseArray()) == ForwardArray()
    @test reverse(ForwardArray()) == ReverseArray()
    @test reverse(ReverseRelation()) == ForwardRelation()
    @test reverse(ForwardRelation()) == ReverseRelation()
    @test reverse(ArrayOrder, Unordered(ForwardRelation())) ==
        Unordered(ReverseRelation())
    @test reverse(ArrayOrder, Ordered(ForwardIndex(), ReverseArray(), ForwardRelation())) ==
        Ordered(ForwardIndex(), ForwardArray(), ReverseRelation())
    @test reverse(IndexOrder, Unordered(ForwardRelation())) ==
        Unordered(ReverseRelation())
    @test reverse(IndexOrder, Ordered(ForwardIndex(), ReverseArray(), ForwardRelation())) ==
        Ordered(ReverseIndex(), ReverseArray(), ReverseRelation())
    @test order(reverse(IndexOrder, Sampled(order=Ordered(ForwardIndex(), ReverseArray(), ForwardRelation()), span=Regular(1)))) ==
        Ordered(ReverseIndex(), ReverseArray(), ReverseRelation())
        Ordered(ForwardIndex(), ForwardArray(), ReverseRelation())
    @test order(reverse(IndexOrder, Sampled(order=Ordered(ForwardIndex(), ReverseArray(), ReverseRelation()), span=Regular(1)))) ==
        Ordered(ReverseIndex(), ReverseArray(), ForwardRelation())
    @test order(reverse(IndexOrder, Categorical(order=Ordered(ForwardIndex(), ReverseArray(), ReverseRelation())))) ==
        Ordered(ReverseIndex(), ReverseArray(), ForwardRelation())
end

@testset "_slicemode" begin
    ind = [10.0, 20.0, 30.0, 40.0, 50.0]

    @testset "Irregular forwards" begin
        mode_ = Sampled(order=Ordered(), span=Irregular((10.0, 60.0)), sampling=Intervals(Start()))
        @test _bounds(_slicemode(mode_, ind, 3), X(ind)) == (30.0, 40.0)
        @test _bounds(_slicemode(mode_, ind, 1:5), X(ind)) == (10.0, 60.0)
        @test _bounds(_slicemode(mode_, ind, 2:3), X(ind)) == (20.0, 40.0)
        mode_ = Sampled(Ordered(), Irregular((0.0, 50.0)), Intervals(End()))
        @test _bounds(_slicemode(mode_, ind, 3), X(ind)) == (20.0, 30.0)
        @test _bounds(_slicemode(mode_, ind, 1:5), X(ind)) == (0.0, 50.0)
        @test _bounds(_slicemode(mode_, ind, 2:3), X(ind)) == (10.0, 30.0)
        mode_ = Sampled(Ordered(), Irregular((5.0, 55.0)), Intervals(Center()))
        @test _bounds(_slicemode(mode_, ind, 3), X(ind)) == (25.0, 35.0)
        @test _bounds(_slicemode(mode_, ind, 1:5), X(ind)) == (5.0, 55.0)
        @test _bounds(_slicemode(mode_, ind, 2:3), X(ind)) == (15.0, 35.0)
    end

    @testset "Irregular reverse" begin
        revind = [50.0, 40.0, 30.0, 20.0, 10.0]
        mode_ = Sampled(order=Ordered(index=ReverseIndex()), span=Irregular(0.0, 50.0),
                        sampling=Intervals(Start()))
        @test _bounds(_slicemode(mode_, revind, 1:5), X(revind)) == (0.0, 50.0)
        @test _bounds(_slicemode(mode_, revind, 1:2), X(revind)) == (30.0, 50.0)
        @test _bounds(_slicemode(mode_, revind, 2:3), X(revind)) == (20.0, 40.0)
        mode_ = Sampled(Ordered(index=ReverseIndex()), Irregular(10.0, 60.0), Intervals(End()))
        @test _bounds(_slicemode(mode_, revind, 1:5), X(revind)) == (10.0, 60.0)
        @test _bounds(_slicemode(mode_, revind, 1:2), X(revind)) == (40.0, 60.0)
        @test _bounds(_slicemode(mode_, revind, 2:3), X(revind)) == (30.0, 50.0)
        mode_ = Sampled(Ordered(index=ReverseIndex()), Irregular(0.5, 55.0), Intervals(Center()))
        @test _bounds(_slicemode(mode_, revind, 1:5), X(revind)) == (0.5, 55.0)
        @test _bounds(_slicemode(mode_, revind, 1:2), X(revind)) == (35.0, 55.0)
        @test _bounds(_slicemode(mode_, revind, 2:3), X(revind)) == (25.0, 45.0)
    end

    @testset "Irregular with no bounds" begin
        mode = Sampled(span=Irregular(nothing, nothing), sampling=Intervals(Start()))
        mode = Sampled(Ordered(), Irregular(nothing, nothing), Intervals(Start()))
        @test _bounds(_slicemode(mode, ind, 3), X()) == (30.0, 40.0)
        @test _bounds(_slicemode(mode, ind, 2:4), X()) == (20.0, 50.0)
        # TODO should this be built into `identify` to at least get one bound?
        @test _bounds(_slicemode(mode, ind, 1:5), X()) == (10.0, nothing)
        mode = Sampled(span=Irregular(nothing, nothing), sampling=Intervals(End()))
        mode = Sampled(Ordered(), Irregular(nothing, nothing), Intervals(End()))
        @test _bounds(_slicemode(mode, ind, 3), X()) == (20.0, 30.0)
        @test _bounds(_slicemode(mode, ind, 2:4), X()) == (10.0, 40.0)
        @test _bounds(_slicemode(mode, ind, 1:5), X()) == (nothing, 50.0)
        mode = Sampled(Ordered(), Irregular(nothing, nothing), Intervals(Center()))
        @test _bounds(_slicemode(mode, ind, 3), X()) == (25.0, 35.0)
        @test _bounds(_slicemode(mode, ind, 2:4), X()) == (15.0, 45.0)
        @test _bounds(_slicemode(mode, ind, 1:5), X()) == (nothing, nothing)
    end

    @testset "Regular is unchanged" begin
        mode = Sampled(span=Regular(1.0), sampling=Intervals(Start()))
        mode = Sampled(Ordered(), Regular(1.0), Intervals(Start()))
        @test _slicemode(mode, ind, 2:3) === mode
    end

    @testset "Points is unchanged" begin
        mode = Sampled(span=Regular(1.0), sampling=Points())
        mode = Sampled(Ordered(), Regular(1.0), Points())
        @test _slicemode(mode, ind, 2:3) === mode
    end

end

@testset "bounds" begin

    @testset "Intervals" begin
        @testset "Regular bounds are calculated from interval type and span value" begin
            @testset "forward ind" begin
                ind = 10.0:10.0:50.0
                dim = X(ind; mode=Sampled(order=Ordered(), sampling=Intervals(Start()), span=Regular(10.0)))
                @test bounds(dim) == (10.0, 60.0)
                dim = X(ind; mode=Sampled(order=Ordered(), sampling=Intervals(End()), span=Regular(10.0)))
                @test bounds(dim) == (0.0, 50.0)
                dim = X(ind; mode=Sampled(Ordered(), Regular(10.0), Intervals(Start())))
                @test bounds(dim) == (10.0, 60.0)
                dim = X(ind; mode=Sampled(Ordered(), Regular(10.0), Intervals(End())))
                @test bounds(dim) == (0.0, 50.0)
                dim = X(ind; mode=Sampled(Ordered(), Regular(10.0), Intervals(Center())))
                @test bounds(dim) == (5.0, 55.0)
            end
            @testset "reverse ind" begin
                revind = [10.0, 9.0, 8.0, 7.0, 6.0]
                dim = X(revind; mode=Sampled(; order=Ordered(ReverseIndex(),ForwardArray(),ForwardRelation()),
                                                   sampling=Intervals(Start()), span=Regular(-1.0)))
                dim = X(revind; mode=Sampled(Ordered(ReverseIndex(),ForwardArray(),ForwardRelation()), Regular(-1.0), Intervals(Start())))
                @test bounds(dim) == (6.0, 11.0)
                dim = X(revind; mode=Sampled(; order=Ordered(ReverseIndex(),ForwardArray(),ForwardRelation()),
                                                   sampling=Intervals(End()), span=Regular(-1.0)))
                dim = X(revind; mode=Sampled(Ordered(ReverseIndex(),ForwardArray(),ForwardRelation()), Regular(-1.0), Intervals(End())))
                @test bounds(dim) == (5.0, 10.0)
                dim = X(revind; mode=Sampled(; order=Ordered(ReverseIndex(),ForwardArray(),ForwardRelation()),
                                                   sampling=Intervals(Center()), span=Regular(-1.0)))
                dim = X(revind; mode=Sampled(Ordered(ReverseIndex(),ForwardArray(),ForwardRelation()), Regular(-1.0), Intervals(Center())))
                @test bounds(dim) == (5.5, 10.5)
            end
        end
        @testset "Irregular bounds are whatever is stored in span" begin
            ind = 10.0:10.0:50.0
            dim = X(ind; mode=Sampled(Ordered(), Irregular(0.0, 50000.0), Intervals(Start())))
            @test bounds(dim) == (0.0, 50000.0)
            @test bounds(DimensionalData._slicedims(getindex, dim, 2:3)[1][1]) == (20.0, 40.0)
        end
        @testset "Explicit bounds are is stored in span matrix" begin
            ind = 10.0:10.0:50.0
            bnds = vcat(ind', (ind .+ 10)')
            dim = X(ind; mode=Sampled(Ordered(), Explicit(bnds), Intervals(Start())))
            @test bounds(dim) == (10.0, 60.0)
            @test bounds(DimensionalData._slicedims(getindex, dim, 2:3)[1][1]) == (20.0, 40.0)
        end
    end

    @testset "Points" begin
        ind = 10:15
        dim = X(ind; mode=Sampled(order=Ordered(), sampling=Points()))
        @test bounds(dim) == (10, 15)
        ind = 15:-1:10
        dim = X(ind; mode=Sampled(order=Ordered(index=ReverseIndex()), sampling=Points()))
        last(dim), first(dim)
        @test bounds(dim) == (10, 15)
        dim = X(ind; mode=Sampled(order=Unordered(), sampling=Points()))
        @test bounds(dim) == (nothing, nothing)
    end

    @testset "Categorical" begin
        ind = [:a, :b, :c, :d]
        dim = X(ind; mode=Categorical(; order=Ordered()))
        @test order(dim) == Ordered()
        @test arrayorder(dim) == ForwardArray()
        @test indexorder(dim) == ForwardIndex()
        @test relation(dim) == ForwardRelation()
        @test_throws ErrorException step(dim)
        @test span(dim) == DimensionalData.NoSpan()
        @test sampling(dim) == DimensionalData.NoSampling()
        @test dims(mode(dim)) === nothing
        @test locus(dim) == Center()
        @test bounds(dim) == (:a, :d)
        dim = X(ind; mode=Categorical(; order=Ordered(; index=ReverseIndex())))
        @test bounds(dim) == (:d, :a)
        @test arrayorder(dim) == ForwardArray()
        @test indexorder(dim) == ReverseIndex()
        @test relation(dim) == ForwardRelation()
        dim = X(ind; mode=Categorical(; order=Unordered()))
        @test bounds(dim) == (nothing, nothing)
        @test order(dim) == Unordered()
    end

end


text:
# cross ratio for 4 numbers
function crossratio(a,b,c,d)
    d_ab = norm(a - b)
    d_bd = norm(b - d)
    d_ac = norm(a - c)
    d_cd = norm(c - d)
    (d_ab/d_bd)/(d_ac/d_cd)
end
export crossratio

# modulo that starts from one
# imod(1,3) = 1
# imod(2,3) = 2
# imod(3,3) = 3
# imod(4,3) = 1
function imod(i, n)
  return (i-1)%n + 1
end
export imod

#TODO handle zero deviation
# l1 = [(x1,y1), (x2,y2)]
# l2 = [(x3,y3), (x4,y4)]
function intersection_point(l1, l2)
  x1 = l1[1][1]
  y1 = l1[1][2]

  x2 = l1[2][1]
  y2 = l1[2][2]

  x3 = l2[1][1]
  y3 = l2[1][2]

  x4 = l2[2][1]
  y4 = l2[2][2]

  denominator = (x1-x2)*(y3-y4)-(y1-y2)*(x3-x4)
  x = (x1*y2-y1*x2)*(x3-x4)-(x1-x2)*(x3*y4-y3*x4)
  y = (x1*y2-y1*x2)*(y3-y4)-(y1-y2)*(x3*y4-y3*x4)
  return [x,y]/denominator
end
export intersection_point

"""
    calcualte_rigid(source_pts,target_pts)

find (R,T) rigid transformation matrix and translation vector
such that target_pts ≈ R * source_pts + T
in the least-square sense
"""
function rigid_trans(src_pts,dst_pts)
  to_array(vv) = [vv[i][j] for i=1:length(vv) , j=1:length(vv[1])]
  n = length(src_pts)
  dim = size(src_pts[1])[1]
  dst_bar = mean(dst_pts)
  src_bar = mean(src_pts)
  dst_centered = dst_pts .- [dst_bar]
  src_centered = src_pts .- [src_bar]

  C = to_array(src_centered)' * to_array(dst_centered) / n
  U,s,V = svd(C)
  # Rotation R in least squares sense:
  Rr = (U * diagm(vcat(ones(dim-1), det(U*V.'))) * V.' ).'
  Tr = -Rr*src_bar + dst_bar

  Pr = PM(Rr,Tr)
  reproj_error = mean(norm.([Rr] .* src_pts .+ [Tr] - dst_pts))
  (Rr,Tr,reproj_error)
end
export rigid_trans


text:
using Flux
using Flux: @epochs

mutable struct NeuralSurrogate{X,Y,M,L,O,P,N,A,U} <: AbstractSurrogate
    x::X
    y::Y
    model::M
    loss::L
    opt::O
    ps::P
    n_echos::N
    lb::A
    ub::U
 end

 """
 NeuralSurrogate(x,y,lb,ub,model,loss,opt,n_echos)

 - model: Flux layers
 - loss: loss function
 - opt: optimization function

 """
 function NeuralSurrogate(x,y,lb,ub; model = Chain(Dense(length(x[1]),1), first), loss = (x,y) -> Flux.mse(model(x), y),opt = Descent(0.01),n_echos::Int = 1)
     X = vec.(collect.(x))
     data = zip(X, y)
     ps = Flux.params(model)
     @epochs n_echos Flux.train!(loss, ps, data, opt)
     return NeuralSurrogate(x,y,model,loss,opt,ps,n_echos,lb,ub)
 end

 function (my_neural::NeuralSurrogate)(val)
     v = [val...]
     out = my_neural.model(v)
     if length(out) == 1
         return out[1]
     else
         return out
     end
 end

function add_point!(my_n::NeuralSurrogate, x_new, y_new)
    if eltype(x_new) == eltype(my_n.x)
        append!(my_n.x, x_new)
        append!(my_n.y, y_new)
    else
        push!(my_n.x, x_new)
        push!(my_n.y, y_new)
    end
    X = vec.(collect.(my_n.x))
    data = zip(X, my_n.y)
    @epochs my_n.n_echos Flux.train!(my_n.loss, my_n.ps, data, my_n.opt)
    nothing
end


text:
abstract type AbstractPartition{Ti} end

struct SplitPartition{Ti} <: AbstractPartition{Ti}
    K::Int
    spl::Vector{Ti}
end

Base.length(Π::SplitPartition) = Π.K

struct DomainPartition{Ti} <: AbstractPartition{Ti}
    K::Int
    prm::Vector{Ti}
    spl::Vector{Ti}
end

Base.length(Π::DomainPartition) = Π.K

struct MapPartition{Ti} <: AbstractPartition{Ti}
    K::Int
    asg::Vector{Ti}
end

Base.length(Π::MapPartition) = Π.K

function Base.:(==)(Π::SplitPartition, Π′::SplitPartition)
    return Π.spl == Π′.spl
end

function Base.:(==)(Π::DomainPartition, Π′::DomainPartition)
    return Π.prm == Π′.prm && Π.spl == Π′.spl
end

function Base.:(==)(Π::MapPartition, Π′::MapPartition)
    return Π.K == Π′.K && Π.asg == Π′.asg
end

function Base.convert(T::Type{<:DomainPartition}, Π::SplitPartition{Ti}) where {Ti}
    prm = collect(1:(Π.spl[end]-1))
    return T(Π.K, prm, Π.spl)
end

function Base.convert(T::Type{<:DomainPartition}, Π::MapPartition{Ti}) where {Ti}
    n = length(Π.asg)
    spl = zeros(Ti, Π.K + 1)
    prm = undefs(Ti, n)
    for j = 1:n
        k = Π.asg[j]
        spl[k + 1] += 1
    end
    q = 1
    for k = 1:(Π.K + 1)
        (spl[k], q) = (q, q + spl[k])
    end
    for j = 1:n
        k = Π.asg[j]
        prm[spl[k + 1]] = j
        spl[k + 1] += 1
    end
    return T(Π.K, prm, spl)
end

function Base.convert(T::Type{<:MapPartition}, Π::SplitPartition{Ti}) where {Ti}
    asg = undefs(Ti, Π.spl[end]-1)
    for k = 1:Π.K
        for j = Π.spl[k]:(Π.spl[k + 1] - 1)
            asg[j] = k
        end
    end
    return T(Π.K, asg)
end

function Base.convert(T::Type{<:MapPartition}, Π::DomainPartition{Ti}) where {Ti}
    asg = undefs(Ti, Π.spl[end]-1)
    for k = 1:Π.K
        for q = Π.spl[k]:(Π.spl[k + 1] - 1)
            j = Π.prm[q]
            asg[j] = k
        end
    end
    return T(Π.K, asg)
end


text:
# function invC(x::Vector{Float64}, Q::AbstractSparseMatrix{T})
#     v = x
#     for j in 1:Q.n
#         tmpx = 0.0
#         for z in Q.colptr[j]:(Q.colptr[j+1]-1)
#             i = Q.rowval[z]
#             if i == j
#                 tmpx += v[j]
#                 v[j] = tmpx / (-Q.nzval[z])
#                 break
#             else
#                 tmpx -= (-Q.nzval[z]) * v[i]
#             end
#         end
#     end
#     v
# end

# function ex(x::Vector{Float64}, Q::SparseMatrixCSC;
#     maxiter=5000, tol=1.0e-16)::Vector{Float64}
#     b = zeros(Q.n)
#     v = x
#     y = invC(v, Q)
#     while maximum(v) > tol
#         gsstep!(v, Q, b)
#         y += v
#     end
#     y
# end

# function ex2(Q::SparseMatrixCSC, b::Vector{Float64};
#     x0::Vector{Float64}=fill(1/Q.n, Q.n),
#     maxiter=5000, steps=50, rtol=1.0e-6)::Vector{Float64}
#     x = x0
#     iter = 0
#     conv = false
#     rerror = 0.0
#     while true
#         prevx = x
#         for i in 1:steps
#             gsstep!(x, Q, b)
#         end
#         rerror = maximum(abs.((x - prevx) ./ x))
#         iter += steps
#         if rerror < rtol
#             conv = true
#             break
#         end
#         if iter >= maxiter
#             break
#         end
#     end
#     @printf "convergence   : %s\n" conv ? "true" : "false"
#     @printf "iteration     : %d / %d\n" iter maxiter
#     @printf "relative error: %e < %e\n" rerror rtol
#     x
# end

# function test2(Q, x)
#     C = triu(Q, 0)
#     D = tril(Q, -1)
#     y = zeros(size(x))
#     x = (-C)' \ x
#     y += x
#     for i in 1:100
#         x = D' * x
#         x = (-C)' \ x
#         y += x
#     end
#     y
# end

# function test(Q::SparseMatrixCSC;
#     x0::Vector{Float64}=fill(1/Q.n, Q.n),
#     maxiter=5000, steps=50, rtol=1.0e-6)::Vector{Float64}
#     x = x0
#     iter = 0
#     conv = false
#     rerror = 0.0
#     while true
#         prevx = x
#         for i in 1:steps
#             b = x
#             gsstep!(x, Q, b)
#             x /= sum(x)
#         end
#         rerror = maximum(abs.((x - prevx) ./ x))
#         iter += steps
#         if rerror < rtol
#             conv = true
#             break
#         end
#         if iter >= maxiter
#             break
#         end
#     end
#     @printf "convergence   : %s\n" conv ? "true" : "false"
#     @printf "iteration     : %d / %d\n" iter maxiter
#     @printf "relative error: %e < %e\n" rerror rtol
#     x
# end


text:
export BERL_CGP, process

function BERL_CGP(cfg::Dict, fitness::Function)::Cambrian.Evolution
    eval = evo -> Cambrian.fitness_evaluate!(evo; fitness=fitness)

    e = Cambrian.Evolution(
        CGPInd,
        cfg;
        id=cfg["id"],
        populate = cgp_populate!,
        evaluate = eval,
    )
    e
end

function cgp_populate!(evo::Cambrian.Evolution)
    mutation = i::CGPInd -> goldman_mutate(evo.cfg, i)
    Cambrian.oneplus_populate!(evo; mutation = mutation, reset_expert = true)
end

function process(i::CGPInd, X)
    CartesianGeneticProgramming.process(i, X)
end

algorithms["CGP"]=BERL_CGP


text:
import OpenCL
using OpenCL: cl

const ARRAY_SIZE = 64
sum = Vector{Float32}(2)
data = Vector{Float32}(ARRAY_SIZE)
for i in 1:ARRAY_SIZE
   data[i] = i*1.0
end

device, ctx, queue = cl.create_compute_context()

# create opencl buffer objects
# copies to the device initiated when the kernel function is called
global_size = 8
local_size = 4
num_groups = global_size/local_size
inputBuffer = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=data)
sumBuffer = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=sum)

addNumbersSource = open(readall, "addNumbers.cl")
# build the program and construct a kernel object
p = cl.Program(ctx, source=addNumbersSource) |> cl.build!
#p = cl.build!(p, options="-w")
const KERNEL_FUNC = "add_numbers"
addNumbersKernel = cl.Kernel(p, KERNEL_FUNC)

LocalMem(::Type{T}, len::Integer) where {T} = begin
    @assert len > 0
    nbytes = sizeof(T) * len
    return LocalMem{T}(convert(Csize_t, nbytes))
end

localBuffer = cl.LocalMem(Float32, local_size)
# call the kernel object with global size set to the size our arrays
addNumbersKernel[queue, (global_size,), (local_size,)](inputBuffer, localBuffer, sumBuffer)

# perform a blocking read of the result from the device
r = cl.read(queue, sumBuffer)

total = 0.0;
for j in 1:length(r)
     total += r[j]
end
actualTotal = (ARRAY_SIZE + 1) * ARRAY_SIZE / 2

println("Total:", total, " Expected:",actualTotal)
println("Norm:", norm(actualTotal - total))
# check to see if our result is what we expect!
if isapprox(norm(actualTotal - total), zero(Float32))
    info("Success!")
else
    error("Norm should be 0.0f")
end


text:
using HDF5Utils
using HDF5
using Glob
using Test

cd(mktempdir())

dict = Dict(
    "str" => "abcdef",
    "mlstr" => s10"北京欢迎你",
    "nt" => (α = 1, β = 2, γ = s3"β"),
    "longstr" => string(rand(100)),
    "arr[uint8]" => zeros(UInt8, 7),
    "arr[mlstr]" => MLString{16}["α" "β"; "γ" "δ"] |> vec,
    "arr[nt]" => repeat([(a = 1, b = s4"α", c = 1f0)], 2),
    "group" => Dict(
        "arr[float32]" => zeros(Float32, 10),
        "group" => Dict("arr[float64]" => ones(10)),
        )
    )
h5save("test.h5", dict)
@test h5open(read, "test.h5") == dict
@test h5open(tryreadmmap, "test.h5") == dict
@test h5load("test.h5") == dict

mutable struct Data
    x::Array{Float32, 2}
    y::Array{Float64, 3}
    z::Array{Float64, 3}
    a::Dict{String, Int}
    b::String
end
x = rand(Float32, 20, 30)
y = rand(20, 30, 40)
z = rand(20, 30, 40)
data = Data(x, y, z, Dict("a" => 1), "abcd")
h5save("data.h5", data)
data′ = h5load("data.h5", Data)
for s in fieldnames(Data)
    @test getfield(data, s) == getfield(data′, s)
end
data′ = 0

for v in (false, true)
    file = "vconcat_$v.h5"
    h5concat("vconcat_$v.h5", ["data.h5", "data.h5"], dims = -2, virtual = v)
    h5concat!("vconcat_$v.h5", "q", r"y|z"; dims = 1, virtual = v)
    @test h5loadv("vconcat_$v.h5", "y") ≈ repeat(y, outer = (1, 2, 1))
    @test h5loadv("vconcat_$v.h5", "q") ≈ cat([repeat(a, outer = (1, 2, 1)) for a in (y, z)]..., dims = 1)
    srcs_list = [["data.h5", "data.h5"], ["data.h5"], ["data.h5", "data.h5", "data.h5"]]
    h5concat("vconcat2d_$v.h5", srcs_list, dims = (-2, -1), virtual = v)
    @test h5loadv("vconcat2d_$v.h5", "x") ≈ [x x zero(x); x zero(x) zero(x); x x x]
end

h5open("diskarray.h5", "w") do fid
    x = rand(10, 9, 8, 7)
    fid["x", "compress", 3] = x
    x′ = HDF5DiskArray(fid["x"])
    @test Array(x′) == x
    @test reshape(x′, 10, 8, 9, 7) == reshape(x, 10, 8, 9, 7)
    @test reshape(x′, 7, 8, 9, 10) == reshape(x, 7, 8, 9, 10)
    @test x′[:, :, 1, 1] == x[:, :, 1, 1]
    @test x′[:, :, 1:2:5, 2:2:6] == x[:, :, 1:2:5, 2:2:6]
    @test x′[:, :, 1:2:5, [1, 2, 5]] == x[:, :, 1:2:5, [1, 2, 5]]
end

GC.gc(true)
h5open("vds.h5", "w") do fid
    layout = VirtualLayout((10, sum(1:5)), Float64)
    for i in 1:5
        h5open("vds_$i.h5", "w") do fid′
            fid′["x"] = rand(10, i)
        end
        off = sum(1:i-1)
        layout[:,  (off + 1):(off + i)] = VirtualSource("vds_$i.h5", "x")[:, :]
    end
    d_create_virtual(fid, "x", layout)
end
h5concat("vdscat.h5", ["vds_$i.h5" for i in 1:5], dims = -1)
@test h5loadv("vdscat.h5", "x") ≈ hcat([h5loadv(h5, "x") for h5 in glob("vds_*.h5")]...)
@test h5loadv("vds.h5", "x") ≈ hcat([h5loadv(h5, "x") for h5 in glob("vds_*.h5")]...)

function sumloop(x)
    s = 0.0
    @inbounds for i in eachindex(x)
        s += x[i]
    end
    return s
end

for d in 1:3
    for n in (10, 100, 200, 500, 10^7)
        GC.gc(true)
        exp(d * log(n)) * 8 > 1024^3 && continue
        h5open("disk_$(d)_$n.h5", "w") do fid
            chunk = ntuple(i -> min(n, 100), d)
            x = rand(fill(n, d)...)
            fid["x", "chunk", chunk] = x
        end
        x = h5load("disk_$(d)_$n.h5", "x")
        sum(read(x.ds)), sumloop(x)
        t_cache = @elapsed s_cache = sumloop(x)
        t_read = @elapsed s_read = sum(read(x.ds))
        @show n, d, t_cache / t_read
        @test s_read ≈ s_cache
    end
end

h5open("attrs.h5", "w") do fid
    fid["x"] = [1]
    fid["y"] = "aaaaaaaa"
    attrs(fid["x"])["name"] = "abc"
end

h5open("attrs.h5", "r") do fid
    @test read(fid["y"]) == "aaaaaaaa"
    a = attrs(fid["x"])["name"]
    @test read(a, String) == "abc"
    @test read(a) == "abc"
    @test read(a) isa MLString
end

!Sys.iswindows() && foreach(rm, glob("*.h5"))


text:
#potrf
for (bname, fname,elty) in ((:cusolverDnSpotrf_bufferSize, :cusolverDnSpotrf, :Float32),
                            (:cusolverDnDpotrf_bufferSize, :cusolverDnDpotrf, :Float64),
                            (:cusolverDnCpotrf_bufferSize, :cusolverDnCpotrf, :Complex64),
                            (:cusolverDnZpotrf_bufferSize, :cusolverDnZpotrf, :Complex128))
    @eval begin
        function potrf!(uplo::BlasChar,
                        A::CudaMatrix{$elty})
            cuuplo = cublasfill(uplo)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("Cholesky factorization is only possible for square matrices!"))
            end
            lda     = max(1,stride(A,2))
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, cublasFillMode_t, Cint,
                               Ptr{$elty}, Cint, Ptr{Cint}),
                              cusolverDnhandle[1], cuuplo, n, A, lda, bufSize))

            buffer  = CudaArray(zeros($elty,bufSize[1]))
            devinfo = CudaArray(zeros(Cint,1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, cublasFillMode_t, Cint,
                               Ptr{$elty}, Cint, Ptr{$elty}, Cint, Ptr{Cint}),
                              cusolverDnhandle[1], cuuplo, n, A, lda, buffer,
                              bufSize[1], devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(-info[1])th parameter is wrong"))
            elseif info[1] > 0
                throw(Base.LinAlg.SingularException(info[1]))
            end
            A
        end
    end
end

#getrf
for (bname, fname,elty) in ((:cusolverDnSgetrf_bufferSize, :cusolverDnSgetrf, :Float32),
                            (:cusolverDnDgetrf_bufferSize, :cusolverDnDgetrf, :Float64),
                            (:cusolverDnCgetrf_bufferSize, :cusolverDnCgetrf, :Complex64),
                            (:cusolverDnZgetrf_bufferSize, :cusolverDnZgetrf, :Complex128))
    @eval begin
        function getrf!(A::CudaMatrix{$elty})
            m,n     = size(A)
            lda = max(1,stride(A,2))
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{$elty}, Cint,
                               Ptr{Cint}), cusolverDnhandle[1], m, n, A, lda,
                              bufSize))

            buffer  = CudaArray(zeros($elty,bufSize[1]))
            devipiv = CudaArray(zeros(Cint,min(m,n)))
            devinfo = CudaArray(zeros(Cint,1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{$elty},
                               Cint, Ptr{$elty}, Ptr{Cint}, Ptr{Cint}),
                              cusolverDnhandle[1], m, n, A, lda, buffer,
                              devipiv, devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            elseif info[1] > 0
                throw(Base.LinAlg.SingularException(info[1]))
            end
            A, devipiv
        end
    end
end

#geqrf
for (bname, fname,elty) in ((:cusolverDnSgeqrf_bufferSize, :cusolverDnSgeqrf, :Float32),
                            (:cusolverDnDgeqrf_bufferSize, :cusolverDnDgeqrf, :Float64),
                            (:cusolverDnCgeqrf_bufferSize, :cusolverDnCgeqrf, :Complex64),
                            (:cusolverDnZgeqrf_bufferSize, :cusolverDnZgeqrf, :Complex128))
    @eval begin
        function geqrf!(A::CudaMatrix{$elty})
            m,n     = size(A)
            lda = max(1,stride(A,2))
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{$elty}, Cint,
                               Ptr{Cint}), cusolverDnhandle[1], m, n, A,
                              lda, bufSize))

            buffer  = CudaArray(zeros($elty, bufSize[1]))
            devtau  = CudaArray(zeros($elty, min(m, n)))
            devinfo = CudaArray(zeros(Cint, 1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{$elty},
                               Cint, Ptr{$elty}, Ptr{$elty}, Cint, Ptr{Cint}),
                              cusolverDnhandle[1], m, n, A, lda, devtau, buffer,
                              bufSize[1], devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            end
            A, devtau
        end
    end
end

#sytrf
for (bname, fname,elty) in ((:cusolverDnSsytrf_bufferSize, :cusolverDnSsytrf, :Float32),
                            (:cusolverDnDsytrf_bufferSize, :cusolverDnDsytrf, :Float64),
                            (:cusolverDnCsytrf_bufferSize, :cusolverDnCsytrf, :Complex64),
                            (:cusolverDnZsytrf_bufferSize, :cusolverDnZsytrf, :Complex128))
    @eval begin
        function sytrf!(uplo::BlasChar,
                        A::CudaMatrix{$elty})
            cuuplo = cublasfill(uplo)
            n      = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("SymTridiagonal matrix must be square!"))
            end
            lda = max(1,stride(A,2))
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Ptr{$elty}, Cint,
                               Ptr{Cint}), cusolverDnhandle[1], n, A, lda,
                              bufSize))

            buffer  = CudaArray(zeros($elty, bufSize[1]))
            devipiv = CudaArray(zeros(Cint, n))
            devinfo = CudaArray(zeros(Cint, 1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, cublasFillMode_t, Cint,
                               Ptr{$elty}, Cint, Ptr{Cint}, Ptr{$elty}, Cint,
                               Ptr{Cint}), cusolverDnhandle[1], cuuplo, n, A,
                              lda, devipiv, buffer, bufSize[1], devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            elseif info[1] > 0
                throw(Base.LinAlg.SingularException(info[1]))
            end
            A, devipiv
        end
    end
end

#potrs
for (fname,elty) in ((:cusolverDnSpotrs, :Float32),
                     (:cusolverDnDpotrs, :Float64),
                     (:cusolverDnCpotrs, :Complex64),
                     (:cusolverDnZpotrs, :Complex128))
    @eval begin
        function potrs!(uplo::BlasChar,
                        A::CudaMatrix{$elty},
                        B::CudaVecOrMat{$elty})
            cuuplo = cublasfill(uplo)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("Cholesky factorization is only possible for square matrices!"))
            end
            if size(B,1) != n
                throw(DimensionMismatch("first dimension of B, $(size(B,1)), must match second dimension of A, $n"))
            end
            nrhs = size(B,2)
            lda  = max(1,stride(A,2))
            ldb  = max(1,stride(B,2))

            devinfo = CudaArray(zeros(Cint,1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, cublasFillMode_t, Cint, Cint,
                               Ptr{$elty}, Cint, Ptr{$elty}, Cint, Ptr{Cint}),
                              cusolverDnhandle[1], cuuplo, n, nrhs, A, lda, B,
                              ldb, devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            end
            B
        end
    end
end

#getrs
for (fname,elty) in ((:cusolverDnSgetrs, :Float32),
                     (:cusolverDnDgetrs, :Float64),
                     (:cusolverDnCgetrs, :Complex64),
                     (:cusolverDnZgetrs, :Complex128))
    @eval begin
        function getrs!(trans::BlasChar,
                        A::CudaMatrix{$elty},
                        ipiv::CudaVector{Cint},
                        B::CudaVecOrMat{$elty})
            cutrans = cublasop(trans)
            n = size(A,1)
            if size(A,2) != n
                throw(DimensionMismatch("LU factored matrix A must be square!"))
            end
            if size(B,1) != n
                throw(DimensionMismatch("first dimension of B, $(size(B,1)), must match second dimension of A, $n"))
            end
            nrhs = size(B,2)
            lda  = max(1,stride(A,2))
            ldb  = max(1,stride(B,2))

            devinfo = CudaArray(zeros(Cint,1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, cublasOperation_t, Cint, Cint,
                               Ptr{$elty}, Cint, Ptr{Cint}, Ptr{$elty}, Cint,
                               Ptr{Cint}), cusolverDnhandle[1], cutrans, n, nrhs,
                              A, lda, ipiv, B, ldb, devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            end
            B
        end
    end
end

#ormqr
for (bname, fname,elty) in ((:cusolverDnSgeqrf_bufferSize, :cusolverDnSormqr, :Float32),
                            (:cusolverDnDgeqrf_bufferSize, :cusolverDnDormqr, :Float64),
                            (:cusolverDnCgeqrf_bufferSize, :cusolverDnCunmqr, :Complex64),
                            (:cusolverDnZgeqrf_bufferSize, :cusolverDnZunmqr, :Complex128))
    @eval begin
        function ormqr!(side::BlasChar,
                        trans::BlasChar,
                        A::CudaMatrix{$elty},
                        devtau::CudaVector{$elty},
                        C::CudaVecOrMat{$elty})
            cutrans = cublasop(trans)
            cuside  = cublasside(side)
            m,k     = size(A)
            lda     = max(1,stride(A,2))
            ldc     = max(1,stride(C,2))
            n       = size(C,2)
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{$elty}, Cint,
                               Ptr{Cint}), cusolverDnhandle[1], m, k, A,
                              lda, bufSize))

            buffer  = CudaArray(zeros($elty, bufSize[1]))
            devinfo = CudaArray(zeros(Cint, 1))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, cublasSideMode_t,
                               cublasOperation_t, Cint, Cint, Cint, Ptr{$elty},
                               Cint, Ptr{$elty}, Ptr{$elty}, Cint, Ptr{$elty},
                               Cint, Ptr{Cint}), cusolverDnhandle[1], cuside,
                              cutrans, m, k, n, A, lda, devtau, C, ldc, buffer,
                              bufSize[1], devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            end
            C
        end
    end
end

#gebrd
for (bname, fname, elty, relty) in ((:cusolverDnSgebrd_bufferSize, :cusolverDnSgebrd, :Float32, :Float32),
                                    (:cusolverDnDgebrd_bufferSize, :cusolverDnDgebrd, :Float64, :Float64),
                                    (:cusolverDnCgebrd_bufferSize, :cusolverDnCgebrd, :Complex64, :Float32),
                                    (:cusolverDnZgebrd_bufferSize, :cusolverDnZgebrd, :Complex128, :Float64))
    @eval begin
        function gebrd!(A::CudaMatrix{$elty})
            m,n     = size(A)
            lda     = max(1,stride(A,2))
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{Cint}),
                             cusolverDnhandle[1], m, n, bufSize))

            buffer  = CudaArray(zeros($elty,bufSize[1]))
            devinfo = CudaArray(zeros(Cint,1))
            D       = CudaArray(zeros($relty,min(m,n)))
            E       = CudaArray(zeros($relty,min(m,n)))
            TAUQ    = CudaArray(zeros($elty,min(m,n)))
            TAUP    = CudaArray(zeros($elty,min(m,n)))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{$elty},
                               Cint, Ptr{$relty}, Ptr{$relty}, Ptr{$elty},
                               Ptr{$elty}, Ptr{$elty}, Cint, Ptr{Cint}),
                              cusolverDnhandle[1], m, n, A, lda, D, E, TAUQ,
                              TAUP, buffer, bufSize[1], devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            end
            A, D, E, TAUQ, TAUP
        end
    end
end

for (bname, fname, elty, relty) in ((:cusolverDnSgesvd_bufferSize, :cusolverDnSgesvd, :Float32, :Float32),
                                    (:cusolverDnDgesvd_bufferSize, :cusolverDnDgesvd, :Float64, :Float64),
                                    (:cusolverDnCgesvd_bufferSize, :cusolverDnCgesvd, :Complex64, :Float32),
                                    (:cusolverDnZgesvd_bufferSize, :cusolverDnZgesvd, :Complex128, :Float64))
    @eval begin
        function gesvd!(jobu::BlasChar,
                        jobvt::BlasChar,
                        A::CudaMatrix{$elty})
            m,n     = size(A)
            lda     = max(1,stride(A,2))
            bufSize = Array(Cint,1)
            statuscheck(ccall(($(string(bname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, Cint, Cint, Ptr{Cint}),
                             cusolverDnhandle[1], m, n, bufSize))

            buffer  = CudaArray(zeros($elty, bufSize[1]))
            rbuffer = CudaArray(zeros($relty, 5 * min(m, n)))
            devinfo = CudaArray(zeros(Cint, 1))
            U       = CudaArray(zeros($elty, m, m))
            ldu     = max(1,stride(U,2))
            S       = CudaArray(zeros($relty, min(m, n)))
            Vt      = CudaArray(zeros($elty, n, n))
            ldvt    = max(1,stride(Vt,2))
            statuscheck(ccall(($(string(fname)),libcusolver), cusolverStatus_t,
                              (cusolverDnHandle_t, BlasChar, BlasChar, Cint,
                               Cint, Ptr{$elty}, Cint, Ptr{$relty}, Ptr{$elty},
                               Cint, Ptr{$elty}, Cint, Ptr{$elty}, Cint,
                               Ptr{$relty}, Ptr{Cint}), cusolverDnhandle[1],
                              jobu, jobvt, m, n, A, lda, S, U, ldu, Vt, ldvt,
                              buffer, bufSize[1], rbuffer, devinfo))
            info = to_host(devinfo)
            if info[1] < 0
                throw(ArgumentError("The $(info[1])th parameter is wrong"))
            end
            U, S, Vt
        end
    end
end


text:
const NN_hopping = 2.8
const graphene_lattice_constant = 2.46
const sublattice_shift = -1 / √(3) * graphene_lattice_constant
const graphene_d1 =
    [graphene_lattice_constant / 2, graphene_lattice_constant * √(3) / 2]
const graphene_d2 =
    [-graphene_lattice_constant / 2, graphene_lattice_constant * √(3) / 2]
const UC_area =
    graphene_d1[1] * graphene_d2[2] - graphene_d1[2] * graphene_d2[1] |> abs

@data Sublattice begin
    A
    B
end

function Base.:!(s::Sublattice)
    @match s begin
        A => B
        B => A
    end
end

function Base.isless(s1::Sublattice, s2::Sublattice)
    @match (s1, s2) begin
        (A, B) => true
        _ => false
    end
    end

@data Spin begin
    SpinUp
    SpinDown
end

function Base.:!(s::Spin)
    @match s begin
        SpinUp => SpinDown
        SpinDown => SpinUp
    end
end

function Base.isless(s1::Spin, s2::Spin)
    @match (s1, s2) begin
        (SpinDown, SpinUp) => true
        _ => false
    end
end

"""
    GrapheneCoord(u::Int, v::Int, sublattice::Sublattice)

Lattice coordinate of a carbon atom.

Each coordinate contains the sublattice index `A` or `B`, as well as
the integer coefficients of the two basis vectors
``d\\times(\\pm 1 \\hat{x} + \\sqrt{3}\\hat{y}) / 2``
(`u` for `+`, `v` for `-`), with ``d = 2.46``Å as the lattice constant.
"""
struct GrapheneCoord
    u::Int
    v::Int
    sublattice::Sublattice
end

function Base.show(io::IO, a::GrapheneCoord)
    print(io, "|", a.u, ", ", a.v, ", ", a.sublattice, "⟩")
end

function Base.isless(a1::GrapheneCoord, a2::GrapheneCoord)
    if a1.sublattice != a2.sublattice
        return isless(a1.sublattice, a2.sublattice)
    else
        if a1.v == a2.v
            return isless(a1.u, a2.u)
        else
            return isless(a1.v, a2.v)
        end
    end
end

"""
    GrapheneState(coord::GrapheneCoord, spin::Spin)

Quantum state of an electron in graphene, denoted by
``|u, v, L\\rangle\\otimes |\\sigma\\rangle`` in the drivation.

The state is given by the [`GrapheneCoord`](@ref) of the orbital, as well as
the electronic spin, which can take values `SpinUp` and `SpinDown`.
"""
struct GrapheneState
    coord::GrapheneCoord
    spin::Spin
end

function Base.show(io::IO, a::GrapheneState)
    print(io, a.coord, "⊗|", a.spin, "⟩")
end

function Base.isless(a1::GrapheneState, a2::GrapheneState)
    if a1.spin != a2.spin
        return isless(a1.spin, a2.spin)
    else
        return isless(a1.coord, a2.coord)
    end
end

"""
    crystal_to_cartesian(coord::GrapheneCoord)

Convert a [`GrapheneCoord`](@ref) to a cartesian point with lengths in Å, where
`GrapheneCoord(0, 0, A)` is at the origin.
"""
function crystal_to_cartesian(coord::GrapheneCoord)
    u = coord.u
    v = coord.v
    x = graphene_d1[1] * u + graphene_d2[1] * v
    y = graphene_d1[2] * u + graphene_d2[2] * v

    return ((x, y + (coord.sublattice == B) * sublattice_shift))
end

"""
    graphene_neighbors(atom::GrapheneCoord)

Determine the nearest neighbors of an `atom` and return a vector of the
corresponding [`GrapheneCoord`](@ref)'s.
"""
function graphene_neighbors(atom::GrapheneCoord)
    u = atom.u
    v = atom.v
    s = !atom.sublattice
    return [
        GrapheneCoord(u, v, s),
        GrapheneCoord(u + (-1)^(s == A), v, s),
        GrapheneCoord(u, v + (-1)^(s == A), s),
    ]
end

"""
    graphene_multiple_neighbors(atom::GrapheneCoord, n::Int)

Return an array of [`GrapheneCoord`](@ref)'s obtained by iteratively running
[`graphene_neighbors`](@ref) `n` times, applying it only to the
newly-added [`GrapheneCoord`](@ref)'s from the past iteration. The entries in
the result are unique and sorted by their distance from `atom`.

# Arguments
* `atom`: [`GrapheneCoord`](@ref) from which the iteration begins
* `n`: number of iterations of [`graphene_neighbors`](@ref).
"""
function graphene_multiple_neighbors(atom::GrapheneCoord, idx::Int)
    res = [atom]
    atoms = [atom]
    for ii = 1:idx
        atoms =
            filter(x -> x ∉ res, mapreduce(graphene_neighbors, vcat, atoms)) |>
            unique
        res = vcat(res, atoms)
    end
    res = sort(
        res,
        by=a ->
            norm(crystal_to_cartesian(atom) .- crystal_to_cartesian(a)),
    )
    return res
end

## Propagator
# When computing Ω, occasionally the integrand becomes small enough to give NaN
# The integrand helper functions help to catch these instances
@inline function Ω_Integrand(z, u, v, x::Float64)
    W = ((z / NN_hopping)^2 - 1.0) / (4.0 * cos(x)) - cos(x) |> Complex
    res =
        exp(1.0im * (u - v) * x) / cos(x) *
        ((W - √(W - 1) * √(W + 1))^abs.(u + v)) / (√(W - 1) * √(W + 1))
    return (isnan(res) ? 0.0 + 0.0im : res)
end

@inline function Ω(z, u, v)
    return ((quadgk(
        x -> 2 * Ω_Integrand(z, u, v, x) / (8.0 * π * NN_hopping^2),
        0.0,
        π / 3,
        2 * π / 3,
        π,
        atol=1e-16,
        rtol=1e-4,
    ))[1])
end

@inline function Ωp_Integrand(z, u, v, x::Float64)
    W = ((z / NN_hopping)^2 - 1.0) / (4.0 * cos(x)) - cos(x) |> Complex
    res =
        2 *
        exp(1.0im * (u - v) * x) *
        ((W - √(W - 1) * √(W + 1))^abs.(u + v + 1)) / (√(W - 1) * √(W + 1))
    return (isnan(res) ? 0.0 + 0.0im : res)
end

@inline function Ωp(z, u, v)
    return ((quadgk(
        x -> 2 * Ωp_Integrand(z, u, v, x) / (8.0 * π * NN_hopping^2),
        0.0,
        π / 3,
        2 * π / 3,
        π,
        atol=1e-16,
        rtol=1e-4,
    ))[1])
end

@inline function Ωn_Integrand(z, u, v, x::Float64)
    W = ((z / NN_hopping)^2 - 1.0) / (4.0 * cos(x)) - cos(x) |> Complex
    res =
        2 *
        exp(1.0im * (u - v) * x) *
        ((W - √(W - 1) * √(W + 1))^abs.(u + v - 1)) / (√(W - 1) * √(W + 1))
    return (isnan(res) ? 0.0 + 0.0im : res)
end

@inline function Ωn(z, u, v)
    return ((quadgk(
        x -> 2 * Ωn_Integrand(z, u, v, x) / (8.0 * π * NN_hopping^2),
        0.0,
        π / 3,
        2 * π / 3,
        π,
        atol=1e-16,
        rtol=1e-4,
    ))[1])
end

# The propagator function picks out the correct element of the Ξ matrix based
# on the sublattices of the graphene coordinates.
function graphene_propagator(a_l::GrapheneState, a_m::GrapheneState, z)
    if a_l.spin != a_m.spin
        return 0.0
    else
        u = (a_l.coord).u - (a_m.coord).u
        v = (a_l.coord).v - (a_m.coord).v
        @match ((a_l.coord).sublattice, (a_m.coord).sublattice) begin
            (A, A) => z * Ω(z, u, v)
            (B, B) => z * Ω(z, u, v)
            (A, B) => -NN_hopping * (Ω(z, u, v) + Ωp(z, u, v))
            (B, A) => -NN_hopping * (Ω(z, u, v) + Ωn(z, u, v))
        end
    end
end

        function propagator_matrix(z, States::Vector{GrapheneState})
    precomputed = Dict{Tuple{Int,Int,GrapheneQFT.Sublattice,GrapheneQFT.Sublattice,GrapheneQFT.Spin,GrapheneQFT.Spin,},ComplexF64,}()

    len_states = length(States)
    out = zeros(ComplexF64, len_states, len_states)
    for ii = 1:len_states
        @inbounds for jj = ii:len_states
            key_ = (
                (States[ii].coord).u - (States[jj].coord).u,
                (States[ii].coord).v - (States[jj].coord).v,
                (States[ii].coord).sublattice,
                (States[jj].coord).sublattice,
                States[ii].spin,
                States[jj].spin,
            )
            c = get(precomputed, key_, NaN)
            if isnan(c)
                res = graphene_propagator(States[ii], States[jj], z)
                out[ii, jj] = res
                out[jj, ii] = res
                precomputed[key_] = res
            else
                out[ii, jj] = c
                out[jj, ii] = c
            end
        end
    end
    return out
end


text:
export diffusionReg, wdiffusionReg, wTVReg,anisoTVReg, wdiffusionRegNodal,wTVRegNodal,computeRegularizer,smallnessReg,logBarrier,logBarrierSquared

function computeRegularizer(regFun::Function,mc::Vector,mref::Vector,MInv::AbstractMesh,alpha)
	R,dR,d2R = regFun(mc,mref,MInv)
	return alpha*R, alpha*dR, alpha*d2R
end

function computeRegularizer(regFun::Array{Function},mc::Vector,mref::Array,MInv::AbstractMesh,alpha::Array{Float64})
	numReg = length(regFun)
	if size(mref,2)!=numReg;
		error("computeRegularizer: number of regularizer (=$numReg) does not match number of reference models (=$(size(mref,2))).")
	end
	if length(alpha)!=numReg;
		error("computeRegularizer: number of regularizer (=$numReg) does not match number of alphas (=$(length(alpha))).")
	end

	R = 0.0; dR = zeros(length(mc)); d2R = spzeros(length(mc),length(mc))
	for k=1:numReg
		Rt,dRt,d2Rt = regFun[k](mc,mref[:,k],MInv)
		R   += alpha[k]*Rt
		dR  += alpha[k]*dRt
		d2R += alpha[k]*d2Rt
	end
	return R, dR, d2R
end

"""
	Rc,dR,d2R = diffusionReg(m,mref,M,Iact=1.0)

	Compute diffusion regularizer
		0.5*||GRAD*(m-mref)||_V^2

	Input:
		m     - model
		mref  - reference model
		M     - Mesh
		Iact  - projector on active cells

	Output
		Rc    - value of regularizer
		dR    - gradient w.r.t. m
		d2R   - Hessian
"""
function diffusionReg(m::Vector,mref,M::AbstractMesh;Iact=1.0)
	# Rc = .5* || Grad*m ||^2
	dm   = m .- mref
	Div  = getDivergenceMatrix(M)
	Div  = Iact'*Div   # project to the active cells
	V    = getVolume(M)
	Af   = getFaceAverageMatrix(M)
	d2R  = Div * sdiag(Af'*vec(Vector(diag(V)))) * Div'
	dR   = d2R*dm
	Rc   = 0.5*dot(dm,dR)
	return Rc,dR,d2R
end

"""
	Rc,dR,d2R = smallnessReg(m,mref,M,Iact=1.0)

	Compute smallness regularizer (L2 difference to reference model)

		R(m) = 0.5*||m-mref||_V^2

	Input:
		m     - model
		mref  - reference model
		M     - Mesh

	Output
		Rc    - value of regularizer
		dR    - gradient w.r.t. m
		d2R   - Hessian
"""
function smallnessReg(m::Vector,mref,M::AbstractMesh)
	dm   = m .- mref
	dR   = dm
	Rc   = 0.5*dot(dm,dm)
	return Rc,dR,sparse(1.0I,length(m),length(m))
end

function wdiffusionReg(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
   # Rc = a1*\\Dx(m-mref)||^2 + .. + a3*\\Dz(m-mref)||^2 + a4*|| m -mref ||^2

   if (M.dim==3) && (length(C) == 0)
      alpha1 = 1; alpha2 = 1; alpha3 = 1; alpha4 = 1e-4;
      Wt   = sdiag([alpha1*ones(M.nf[1]);alpha2*ones(M.nf[2]);alpha3*ones(M.nf[3])])
   elseif (M.dim==2) && (length(C) == 0)
	  alpha1 = 1; alpha2 = 1;alpha4 = 1e-4;
	  Wt   = sdiag([alpha1*ones(M.nf[1]);alpha2*ones(M.nf[2])])
   elseif (M.dim==3) && (length(C) == 4)
      # C = alphax, alphay, alphaz, alphas
      alpha1 = C[1]; alpha2 = C[2]; alpha3 = C[3]; alpha4 = C[4]
      Wt   = sdiag([alpha1*ones(M.nf[1]);alpha2*ones(M.nf[2]);alpha3*ones(M.nf[3])])
   elseif length(C) == 5
      #Af  = getFaceAverageMatrix(M)
      #wt  = Af'*C[5]
      alpha1 = C[1]; alpha2 = C[2]; alpha3 = C[3]; alpha4 = C[4]
      wt   = C[5]
      Wt   = sdiag([alpha1*ones(M.nf[1]);alpha2*ones(M.nf[2]);alpha3*ones(M.nf[3])]) *
             sdiag(wt)
      #Wt = sdiag(C[5])
   elseif length(C) == sum(M.nf) + 1
      # C contains alphas followed by interface weights.
      alpha4 = C[1]
      Wt = sdiag( C[2:end] )
      #println("inteface weights used.")
   else
      error("bad regparams in wdiffusionReg")
   end

   dm   = m .- mref
   Div  = getDivergenceMatrix(M)

   V    = getVolume(M)
   #Div  = Iact'*(Div*Wt)   # project to the active cells
   Div  = Iact'*(V*Div*Wt)   # project to the active cells

   Af   = getFaceAverageMatrix(M)

   #d2R  = Div * sdiag(Af'*Vector(diag(V))) * Div' + alpha4*Iact'*V*Iact
   d2R  = Div * sdiag(1 ./(Af'*Vector(diag(V)))) * Div' + alpha4*Iact'*V*Iact

   dR   = d2R*dm
   Rc   = 0.5*dot(dm,dR)

   return Rc,dR,d2R
end # function wdiffusionReg

"""
	Rc,dR,d2R = wTVReg(m,mref,M,Iact,C=[])

	Compute weighted total variation regularizer

	Input:
		m     - model
		mref  - reference model
		M     - Mesh
		Iact  - projector on active cells
		C     - anisotropy parameters (default: [1 1 1])
		eps   - conditioning parameter for TV norm (default: 1e-3)

	Output
		Rc    - value of regularizer
		dR    - gradient w.r.t. m
		d2R   - Hessian
"""
function wTVReg(m::Vector,mref,M::AbstractMesh; Iact=1.0, C=[],eps=1e-3)
	# Rc = sqrt(a1*\\Dx(m-mref)||^2 + .. + a3*\\Dz(m-mref)||^2) + a4*|| m -mref ||^2
	if length(C) == 0
		alpha1 = 1; alpha2 = 1; alpha3 = 1;
	elseif length(C) == 1
		alpha1 = 1; alpha2 = 1; alpha3 = 1;
	else
		alpha1 = C[1]; alpha2 = C[2]; alpha3 = C[3];
	end
	dm   = m .- mref
	Div  = getDivergenceMatrix(M)
	if M.dim==3
		Wt   = sdiag([alpha1*ones(M.nf[1]);alpha2*ones(M.nf[2]);alpha3*ones(M.nf[3])])
	elseif M.dim==2
		Wt   = sdiag([alpha1*ones(M.nf[1]);alpha2*ones(M.nf[2])])
	end
	Div  = Iact'*(Div*Wt)   # project to the active cells
	V    = getVolume(M); v = Vector(diag(V))
	Af   = getFaceAverageMatrix(M)

	wTV  = sqrt.(Af*(Div'*dm).^2 .+eps);

	Rc   = dot(v,wTV);
	d2R  = Div*(sdiag(Af'*(v./wTV))*Div')
	dR   = d2R*dm;
	return Rc,dR,d2R
end

"""
	Rc,dR,d2R = anisoTVReg(m,mref,M,Iact,eps)

	Compute anisotropic total variation regularizer

	Input:
		m     - model
		mref  - reference model
		M     - Mesh
		Iact  - projector on active cells
		eps   - conditioning parameter for TV norm (default: 1e-3)

	Output
		Rc    - value of regularizer
		dR    - gradient w.r.t. m
		d2R   - Hessian
"""
function anisoTVReg(m::Vector,mref,M::AbstractMesh; Iact=1.0, eps=1e-3)
	# Rc = sqrt(a1*\\Dx(m-mref)||^2 + .. + a3*\\Dz(m-mref)||^2) + a4*|| m -mref ||^2
	dm   = m .- mref
	Div  = getVolume(M)*getDivergenceMatrix(M)
	Div  = Iact'*(Div)   # project to the active cells
	v    = ones(size(Div,2))

	wTV  = sqrt.((Div'*dm).^2 .+eps);
	Rc   = dot(v,wTV);
	d2R  = Div*(sdiag(v./wTV)*Div')
	dR   = d2R*dm;
	return Rc,dR,d2R
end

"""
	Rc,dR,d2R = wdiffusionRegNodal(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])

	Computes weighted diffusion regularizer for nodal model

	Input:
		m     - model
		mref  - reference model
		M     - Mesh
		Iact  - projector on active cells
		C     - optional parameters

	Output
		Rc    - value of regularizer
		dR    - gradient w.r.t. m
		d2R   - Hessian
"""
function wdiffusionRegNodal(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
	dm = m.-mref;
	if isempty(C)
		C = [1,1,1,1,1e-5];
	end
	if M.dim==3
		Wt   = [C[1]*ones(M.ne[1]);C[2]*ones(M.ne[2]);C[3]*ones(M.ne[3])];
	else
		Wt   = [C[1]*ones(M.ne[1]);C[3]*ones(M.ne[2])];
	end

	V    = getVolume(M);
	Af   = getEdgeAverageMatrix(M)
	v    = (Af'*Vector(diag(V)))
	Wt   = sdiag(Wt.*(v.*Wt));

	Av = getNodalAverageMatrix(M);

	Av = Av*Iact;
	mass = Av'*V*Av;

	Grad   = getNodalGradientMatrix(M)*Iact;

	d2R = Grad'*Wt*Grad;
	d2R += C[4]*mass;
	dR  = d2R*dm;
	Rc  = 0.5*dot(dm,dR);

	if isnan(Rc)
		dump(m)
	end
	clear!(M);
   return Rc,dR,d2R
end

function wTVRegNodal(m::Vector, mref::Vector, M::AbstractMesh; Iact=1.0, C=[])
	dm = m.-mref;
	if isempty(C)
		C = [1,1,1,1,1e-5];
	end
	eps = 1e-3;
	if M.dim==3
		Wt   = sdiag([C[1]*ones(M.ne[1]);C[2]*ones(M.ne[2]);C[3]*ones(M.ne[3])]);
	else
		Wt   = sdiag([C[1]*ones(M.ne[1]);C[3]*ones(M.ne[2])]);
	end

	V    = getVolume(M);
	Av = getNodalAverageMatrix(M);
	mass = Iact'*Av'*V*Av*Iact;
	d2R  = C[4]*mass;
	Rc   = 0.5*dot(dm,d2R*dm);

	v      = Vector(diag(V))
	Af     = getEdgeAverageMatrix(M)
	Grad   = Wt*getNodalGradientMatrix(M)*Iact;
	wTV    = sqrt.(Af*(((Grad*dm).^2).+eps));
	Rc     += dot(v,wTV);
	d2R    += Grad'*sdiag(Af'*(v./wTV))*Grad;
	dR     = d2R*dm
	if isnan(Rc)
		dump(m)
	end
	clear!(M);
   return Rc,dR,d2R
end

"""
	Rc,dR,d2R = logBarrier(m::Vector, z::Vector, M::AbstractMesh,low::Vector,high::Vector, epsilon)

	Computes logBarrier regularizer

	R = -log(1 - ((m-high)/epsilon).^2) if high-epsilon <  m < high
					0					if low+epsilon  <= m <= high-epsilon
		-log(1 - ((m-low)/epsilon).^2)	if low          <  m < low+epsilon

	Input:
		m     	- model
		z     	- not being used. Here for compatibility.
		M     	- Mesh. not being used. Here for compatibility.
		low   	- low bound for each coordinate.
		high  	- high bound for each coordinate.
		epsilon - layer width of the barier.

	Output
		g    - value of regularizer
		dg    - gradient w.r.t. m
		d2g   - Hessian (diagonal matrix). Second derivative is not continous.
"""
function logBarrier(m::Vector,z::Vector,M::AbstractMesh, low::Vector, high::Vector ,epsilon = min(0.1*abs(low),0.1*abs(high)))

	if length(findall(m.>=high)) + length(findall(m.<=low))>0
		return Inf, zeros(length(m)), spzeros(length(m),length(m));
	end

	z = copy(m);
	indProj = zeros(length(m));

	low = low + epsilon;
	high = high - epsilon;
	t = z.<low;
	indProj[t]  .= 1.0;
	z[t]  = low[t];
	t = z.>high;
	indProj[t] .= 1.0;
	z[t] = high[t];

	e 	= (epsilon.+1e-10);
	dm  = (m-z)./e;
	f   = 1.0 .- dm.^2; # > 0
	if minimum(f) < 0.0
		error("negative f");
	end
	df  = (-2.0).*indProj.*(dm./e);
	d2f = (-2.0).*(indProj./(e.^2));

	g   = -sum(log.(f));
	dg  = -df./f;
	d2g = dg.^2 - d2f./f  ;
	return g,dg,sdiag(d2g);
end

"""
	Rc,dR,d2R = logBarrierSquared(m::Vector, z::Vector, M::AbstractMesh,low::Vector,high::Vector, epsilon)

	Computes logBarrier regularizer

	R = (log(1 - ((m-high)/epsilon).^2))^2 if high-epsilon <  m < high
					0					   if low+epsilon  <= m <= high-epsilon
		(log(1 - ((m-low)/epsilon).^2))^2  if low          <  m < low+epsilon

	Input:
		m     	- model
		z     	- not being used. Here for compatibility.
		M     	- Mesh. not being used. Here for compatibility.
		low   	- low bound for each coordinate.
		high  	- high bound for each coordinate.
		epsilon - layer width of the barier.

	Output
		g    - value of regularizer
		dg    - gradient w.r.t. m
		d2g   - Gauss Newton Hessian approximation (diagonal matrix). Second derivative approx is continous.
"""
function logBarrierSquared(m::Vector,z::Vector,M::AbstractMesh, low::Vector, high::Vector ,epsilon = min(0.1*abs(low),0.1*abs(high)))
	if length(findall(m.>=high)) + length(findall(m.<=low))>0
		return Inf, zeros(length(m)), spzeros(length(m),length(m));
	end
	z = copy(m);
	indProj = zeros(length(m));
	low = low + epsilon;
	high = high - epsilon;
	t = z.<low;
	indProj[t]  .= 1.0;
	z[t]  = low[t];
	t = z.>high;
	indProj[t] .= 1.0;
	z[t] = high[t];

	e 	= (epsilon.+1e-10);
	dm  = (m-z)./e;
	f   = 1.0 .- dm.^2; # > 0
	if minimum(f) < 0.0
		error("negative f");
	end
	df  = (-2.0).*indProj.*(dm./e);

	r   = log.(f);
	dr  = df./f;
	g   = 0.5*dot(r,r);
	dg  = dr.*r;
	d2g = dr.^2;
	return g,dg,sdiag(d2g);
end

export TikhonovReg
function TikhonovReg(m::Vector,mref::Vector,M::AbstractMesh,sigma::SparseMatrixCSC)
d = (m-mref);
df = sigma*d;
f = 0.5*dot(d,df);
d2f = sigma;
return f,df,d2f;
end


text:
using Distributed
addprocs(2, exeflags="--project")

training_step = RemoteChannel(()->Channel{Int}(1))
put!(training_step, 0)
@everywhere begin
function update_remote_counter(remote_counter::RemoteChannel, update_count::Int)
	old_value=take!(remote_counter)
	put!(remote_counter, old_value + update_count)
	return nothing
end

function uses_ts(training_step)

	while fetch(training_step) < 1000
		training_step_ = fetch(training_step)
		println("From user", training_step_)
	end
	return true
end

function changes_ts(training_step)
	while fetch(training_step) < 1000
		update_remote_counter(training_step, 1)
		training_step_ = fetch(training_step)
		println("From changeer", training_step_)
	end
	return true
end
end

u = @spawnat :any uses_ts(training_step)
c = @spawnat :any changes_ts(training_step)


text:
# Iteration occurs as a nested container

# Series
@inline start(s::AbstractSeries) = start(indices(s)[1])
@inline function next(s::AbstractSeries, i)
    (ind, i) = next(indices(s)[1], i)
    @inbounds return (s[ind], i)
end
@inline done(s::AbstractSeries, i) = done(indices(s)[1], i)

@inline length(s::AbstractSeries) = length(indices(s)[1])
@inline endof(s::AbstractSeries) = last(indices(s)[1])

# Tables iterate rows as series
@inline start(t::AbstractTable) = start(indices(t)[1])
@inline function next(t::AbstractTable, i)
    (ind, i) = next(indices(t)[1], i)
    @inbounds return (t[ind, :], i)
end
@inline done(t::AbstractTable, i) = done(indices(t)[1], i)

@inline endof(t::AbstractTable) = last(indices(t)[1]) # not type stable
@inline endof(t::AbstractTable, i::Int) = last(indices(t)[i]) # not type stable
@inline endof(t::AbstractTable, ::Type{Val{i}}) where {i} = last(indices(t)[i::Int])

@inline length(t::AbstractTable) = length(indices(t)[1])


text:
@testset "HamiltonianProblem Test" begin

    @testset "Scalar case test" begin
        H(dθ, θ, p) = dθ / 2 - 9.8 * cos(θ)
        p0, q0   = rand(2)
        acc      = (v, x, p, t) -> ForwardDiff.derivative(x -> -H(v[1], x, p), x[1])
        vel      = (v, x, p, t) -> ForwardDiff.derivative(v -> H(v, x[1], p), v[1])

        prob1    = HamiltonianProblem(H, p0, q0, (0., 10.))
        prob_1   = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))

        @test test_solve(prob1, prob_1)
    end

    @testset "Static Vector case test" begin
        p0 = @SVector rand(2)
        q0 = @SVector rand(2)
        H(dθ, θ, p) = dθ[1] / 2 + dθ[2] / 2 - 9.8 * cos(θ[1]) - 9.8 * cos(θ[2])
        prob1    = HamiltonianProblem(H, p0, q0, (0., 10.))
        acc      = (v, x, p, t) -> ForwardDiff.gradient(x -> -H(v, x, p), x)
        vel      = (v, x, p, t) -> ForwardDiff.gradient(v -> H(v, x, p), v)
        prob_1   = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))
        @test test_solve(prob1, prob_1)
    end

    @testset "Vector case test" begin
        H(dθ, θ, p) = (dθ / 2 - 9.8 * cos.(θ))[1]
        p0, q0 = [rand(5) for i in 1:2]
        acc = (dv, v, x, p, t) -> ForwardDiff.gradient!(dv, x -> -H(v, x, p), x)
        vel = (dx, v, x, p, t) -> ForwardDiff.gradient!(dx, v -> H(v, x, p), v)
        prob2  = HamiltonianProblem(H, p0, q0, (0., 10.))
        prob_2 = DynamicalODEProblem(acc, vel, p0, q0, (0., 10.))
        @test test_solve(prob2, prob_2)
    end
end


text:
"""
```
corners = imcorner(img; [method])
corners = imcorner(img, threshold, percentile; [method])
```

Performs corner detection using one of the following methods -

    1. harris
    2. shi_tomasi
    3. kitchen_rosenfeld

The parameters of the individual methods are described in their documentation. The
maxima values of the resultant responses are taken as corners. If a threshold is
specified, the values of the responses are thresholded to give the corner pixels.
The threshold is assumed to be a percentile value unless `percentile` is set to false.
"""
function imcorner(img::AbstractArray; method::Function = harris, args...)
    responses = method(img; args...)
    corners = similar(img, Bool)
    fill!(corners, false)
    maxima = map(CartesianIndex, findlocalmaxima(responses))
    for m in maxima corners[m] = true end
    corners
end

function imcorner(img::AbstractArray, threshold; method::Function = harris, args...)
    responses = method(img; args...)
    map(i -> i > threshold, responses)
end

function imcorner(img::AbstractArray, thresholdp::Percentile; method::Function = harris, args...)
    responses = method(img; args...)
    threshold = StatsBase.percentile(vec(responses), thresholdp.p)
    map(i -> i > threshold, responses)
end

"""
```
corners = imcorner_subpixel(img; [method])
         -> Vector{HomogeneousPoint{Float64,3}}
corners = imcorner_subpixel(img, threshold, percentile; [method])
         -> Vector{HomogeneousPoint{Float64,3}}
```

Same as [`imcorner`](@ref), but estimates corners to sub-pixel precision.

Sub-pixel precision is achieved by interpolating the corner response values using
the 4-connected neighbourhood of a maximum response value.
See [`corner2subpixel`](@ref) for more details of the interpolation scheme.

"""
function imcorner_subpixel(img::AbstractArray; method::Function = harris, args...)
    responses = method(img; args...)
    corner_indicator = similar(img, Bool)
    fill!(corner_indicator, false)
    maxima = map(CartesianIndex, findlocalmaxima(responses))
    for m in maxima corner_indicator[m] = true end
    corners = corner2subpixel(responses,corner_indicator)
end

function imcorner_subpixel(img::AbstractArray, threshold; method::Function = harris, args...)
    responses = method(img; args...)
    corner_indicator = map(i -> i > threshold, responses)
    corners = corner2subpixel(responses,corner_indicator)
end

function imcorner_subpixel(img::AbstractArray, thresholdp::Percentile; method::Function = harris, args...)
    responses = method(img; args...)
    threshold = StatsBase.percentile(vec(responses), thresholdp.p)
    corner_indicator = map(i -> i > threshold, responses)
    corners = corner2subpixel(responses,corner_indicator)
end

"""
```
corners = corner2subpixel(responses::AbstractMatrix,corner_indicator::AbstractMatrix{Bool})
        -> Vector{HomogeneousPoint{Float64,3}}
```
Refines integer corner coordinates to sub-pixel precision.

The function takes as input a matrix representing corner responses
and a boolean indicator matrix denoting the integer coordinates of a corner
in the image. The output is a vector of type [`HomogeneousPoint`](@ref)
storing the sub-pixel coordinates of the corners.

The algorithm computes a correction factor which is added to the original
integer coordinates. In particular, a univariate quadratic polynomial is fit
separately to the ``x``-coordinates and ``y``-coordinates of a corner and its immediate
east/west, and north/south neighbours. The fit is achieved using a local
coordinate system for each corner, where the origin of the coordinate system is
a given corner, and its immediate neighbours are assigned coordinates of  minus
one and plus one.

The corner and its two neighbours form a system of three equations. For example,
let  ``x_1 = -1``,  ``x_2 = 0`` and  ``x_3 = 1`` denote the local ``x`` coordinates
of the west, center and east pixels and let the vector ``\\mathbf{b} = [r_1, r_2, r_3]``
denote the corresponding corner response values. With

```math
    \\mathbf{A} =
        \\begin{bmatrix}
            x_1^2 & x_1  & 1  \\\\
            x_2^2 & x_2  & 1 \\\\
            x_3^2 & x_3  & 1 \\\\
        \\end{bmatrix},
```
the coefficients of the quadratic polynomial can be found by solving the
system of equations ``\\mathbf{b} = \\mathbf{A}\\mathbf{x}``.
The result is given by ``x = \\mathbf{A}^{-1}\\mathbf{b}``.

The vertex of the quadratic polynomial yields a sub-pixel estimate of the
true corner position. For example, for a univariate quadratic polynomial
``px^2 + qx + r``, the ``x``-coordinate of the vertex is ``\\frac{-q}{2p}``.
Hence, the refined sub-pixel coordinate is equal to:
 ``c +  \\frac{-q}{2p}``, where ``c`` is the integer coordinate.

!!! note
    Corners on the boundary of the image are not refined to sub-pixel precision.

"""
function corner2subpixel(responses::AbstractMatrix, corner_indicator::AbstractMatrix{Bool})
    row_range, col_range = axes(corner_indicator)
    idxs = findall(!iszero, corner_indicator) # findnz
    row, col = (getindex.(idxs,1),getindex.(idxs,2))
    ncorners = length(row)
    corners = fill(HomogeneousPoint((0.0,0.0,0.0)),ncorners)
    invA = @SMatrix [0.5 -1.0 0.5; -0.5 0.0 0.5; 0.0 1.0 -0.0]
    for k = 1:ncorners
        # Corners on the perimeter of the image will not be interpolated.
        if  (row[k] == first(row_range) || row[k] == last(row_range) ||
             col[k] == first(col_range) || col[k] == last(col_range))
            y = convert(Float64,row[k])
            x = convert(Float64,col[k])
            corners[k] = HomogeneousPoint((x,y,1.0))
        else
            center, north, south, east, west =
                                unsafe_neighbourhood_4(responses,row[k],col[k])
            # Solve for the coefficients of the quadratic equation.
            a, b, c = invA* @SVector [west, center, east]
            p, q, r = invA* @SVector [north, center, south]
            # Solve for the first coordinate of the vertex.
            u = -b/(2.0a)
            v = -q/(2.0p)
            corners[k] = HomogeneousPoint((col[k]+u,row[k]+v,1.0))
        end
    end
    return corners
end

"""
```
unsafe_neighbourhood_4(matrix::AbstractMatrix,r::Int,c::Int)
```

Returns the value of a matrix at given coordinates together with the values
of the north, south, east and west neighbours.

This function does not perform bounds checking. It is up to the user to ensure
that the function is not called with indices that are on the boundary of the
matrix.

"""
function unsafe_neighbourhood_4(matrix::AbstractMatrix,r::Int,c::Int)
    center = matrix[r,c]
    north = matrix[r-1,c]
    south = matrix[r+1,c]
    east = matrix[r,c+1]
    west = matrix[r,c-1]
    return center, north, south, east, west
end

"""
```
harris_response = harris(img; [k], [border], [weights])
```

Performs Harris corner detection. The covariances can be taken using either a mean
weighted filter or a gamma kernel.
"""
function harris(img::AbstractArray; k::Float64 = 0.04, args...)
    cov_xx, cov_xy, cov_yy = gradcovs(img, args...)
    corner = map((xx, yy, xy) -> xx * yy - xy ^ 2 - k * (xx + yy) ^ 2, cov_xx, cov_yy, cov_xy)
    corner
end

"""
```
shi_tomasi_response = shi_tomasi(img; [border], [weights])
```

Performs Shi Tomasi corner detection. The covariances can be taken using either a mean
weighted filter or a gamma kernel.
"""
function shi_tomasi(img::AbstractArray; border::AbstractString = "replicate", args...)
    cov_xx, cov_xy, cov_yy = gradcovs(img, border; args...)
    corner = map((xx, yy, xy) -> ((xx + yy) - (sqrt((xx - yy) ^ 2 + 4 * xy ^ 2))) / 2, cov_xx, cov_yy, cov_xy)
    corner
end

"""
```
kitchen_rosenfeld_response = kitchen_rosenfeld(img; [border])
```

Performs Kitchen Rosenfeld corner detection. The covariances can be taken using either a mean
weighted filter or a gamma kernel.
"""
function kitchen_rosenfeld(img::AbstractArray; border::AbstractString = "replicate")
    meth = KernelFactors.sobel
    (grad_x, grad_y) = imgradients(img, meth, border)
    (grad_xx, grad_xy) = imgradients(grad_x, meth, border)
    (grad_yx, grad_yy) = imgradients(grad_y, meth, border)
    map(kr, grad_x, grad_y, grad_xx, grad_xy, grad_yy)
end

function kr(x::T, y::T, xx::T, xy::T, yy::T) where T<:Real
    num = xx*y*y + yy*x*x - 2*xy*x*y
    denom = x*x + y*y
    ifelse(denom == 0, zero(num)/one(denom), -num/denom)
end

function kr(x::Real, y::Real, xx::Real, xy::Real, yy::Real)
    xp, yp, xxp, xyp, yyp = promote(x, y, xx, xy, yy)
    kr(xp, yp, xxp, xyp, yyp)
end

kr(x::RealLike, y::RealLike, xx::RealLike, xy::RealLike, yy::RealLike) =
    kr(gray(x), gray(y), gray(xx), gray(xy), gray(yy))

function kr(x::AbstractRGB, y::AbstractRGB, xx::AbstractRGB, xy::AbstractRGB, yy::AbstractRGB)
    krrgb = RGB(kr(red(x), red(y), red(xx), red(xy), red(yy)),
                kr(green(x), green(y), green(xx), green(xy), green(yy)),
                kr(blue(x),  blue(y),  blue(xx),  blue(xy),  blue(yy)))
    gray(convert(Gray, krrgb))
end

"""
    fastcorners(img, n, threshold) -> corners

Performs FAST Corner Detection. `n` is the number of contiguous pixels
which need to be greater (lesser) than intensity + threshold (intensity - threshold)
for a pixel to be marked as a corner. The default value for n is 12.
"""
function fastcorners(img::AbstractArray{T}, n::Int = 12, threshold::Float64 = 0.15) where T
    img_padded = padarray(img, Fill(0, (3,3)))
    corner = falses(size(img))
    R = CartesianIndices(size(img))
    idx = map(CartesianIndex{2}, [(0, 3), (1, 3), (2, 2), (3, 1), (3, 0), (3, -1), (2, -2), (1, -3),
            (0, -3), (-1, -3), (-2, -2), (-3, -1), (-3, 0), (-3, 1), (-2, 2), (-1, 3)])

    idxidx = [1, 5, 9, 13]
    for I in R
        bright_threshold = img_padded[I] + threshold
        dark_threshold = img_padded[I] - threshold
        if n >= 12
            sum_bright = 0
            sum_dark = 0
            for k in idxidx
                pixel = img_padded[I + idx[k]]
                if pixel > bright_threshold
                    sum_bright += 1
                elseif pixel < dark_threshold
                    sum_dark += 1
                end
            end
            if sum_bright < 3 && sum_dark < 3
                continue
            end
        end
        consecutive_bright = 0
        consecutive_dark = 0

        for i in 1:15 + n
            k = mod1(i, 16)
            pixel = img_padded[I + idx[k]]
            if pixel > bright_threshold
                consecutive_dark = 0
                consecutive_bright += 1
            elseif pixel < dark_threshold
                consecutive_bright = 0
                consecutive_dark += 1
            end

            if consecutive_dark == n || consecutive_bright == n
                corner[I] = true
                break
            end
        end
    end
    corner
end

function gradcovs(img::AbstractArray, border::AbstractString = "replicate"; weights::Function = meancovs, args...)
    (grad_x, grad_y) = imgradients(img, KernelFactors.sobel, border)

    cov_xx = dotc.(grad_x, grad_x)
    cov_xy = dotc.(grad_x, grad_y)
    cov_yy = dotc.(grad_y, grad_y)

    weights(cov_xx, cov_xy, cov_yy, args...)
end

function meancovs(cov_xx, cov_xy, cov_yy, blockSize::Int = 3)

    box_filter_kernel = centered((1 / (blockSize * blockSize)) * ones(blockSize, blockSize))

    filt_cov_xx = imfilter(cov_xx, box_filter_kernel)
    filt_cov_xy = imfilter(cov_xy, box_filter_kernel)
    filt_cov_yy = imfilter(cov_yy, box_filter_kernel)

    filt_cov_xx, filt_cov_xy, filt_cov_yy
end

function gammacovs(cov_xx, cov_xy, cov_yy, gamma::Float64 = 1.4)
    kernel = KernelFactors.gaussian((gamma, gamma))

    filt_cov_xx = imfilter(cov_xx, kernel)
    filt_cov_xy = imfilter(cov_xy, kernel)
    filt_cov_yy = imfilter(cov_yy, kernel)

    filt_cov_xx, filt_cov_xy, filt_cov_yy
end


text:
using Printf

function file_io_examples()
    println("--------------------------------------")
    println("File I/O examples")

    println()

    # Open a file for writing
    open("random.txt", "w") do file
        write(file, "Here is some random text\nIt is great\n")
    end

    # Open a file for reading
    open("random.txt") do file
        # Read whole file into a string
        data = read(file, String)
        println(data)
    end

    open("random.txt") do file
        # Read each line 1 at a time
        for line in eachline(file)
            println(line)
        end
    end
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump