text:
using FFTW
# load data
data_path = "testData/small"
TEs = [4, 8, 12]

mag_nii = readmag("$data_path/Mag.nii")
hdr = mag_nii.header
phase_nii = readphase("$data_path/Phase.nii")

data = Data(mag_nii, phase_nii, hdr, TEs)

# echo formatting test for 2d row vector
Data(mag_nii, phase_nii, hdr, [4 8 12])
print("real data tests")
# default test
swi = calculateSWI(data, Options(phase_unwrap=:romeo))
mip = createMIP(swi)
mIP = createIntensityProjection(swi, minimum)
MIP = createIntensityProjection(swi, maximum)
mean(I) = sum(I) / length(I)
meanIP = createIntensityProjection(swi, mean)

# single-echo
se_data = Data(mag_nii[:,:,:,1], phase_nii[:,:,:,1], hdr)
swi = calculateSWI(data)
mip = createMIP(swi)

options = [
    Options()
    Options(phase_hp_σ=[6,6,3])
    Options(phase_unwrap=:romeo)
    Options(phase_unwrap=:laplacianslice)
    Options(phase_scaling_type=:positive)
    Options(phase_scaling_type=:negativetanh)
    Options(phase_scaling_type=:negative)
    Options(phase_scaling_type=:triangular)
    Options(phase_scaling_strength=6)
    Options(phase_scaling_strength=200)
    Options(mag_combine=:CNR=>(:wm, :gm))
    Options(mag_combine=:CNR=>(:wm, :ven))
    Options(mag_combine=:CNR=>(:wm, :csf))
    Options(mag_combine=:CNR=>(:ms, :iron_ring))
    Options(mag_combine=:average)
    Options(mag_combine=:SE=>10)
    Options(mag_combine=:SE=>5)
    Options(mag_combine=(:closest=>5)) # 1
    Options(mag_combine=(:echo=>2)) # 2
    Options(mag_combine=:last) # 3
    Options(mag_softplus=false)
    Options(mag_sens=[1])
]
# TODO
wrong_options = [
    Options(;phase_scaling_strength=-5)
    Options(;phase_scaling_strength=0)
    Options(;mag_combine=:SE=>0)
    Options(;mag_combine=:SE=>500)
]

s = []
for o in options
    push!(s, calculateSWI(data, o)[:,:,20])
end

# all results should be different
for i in 1:length(s), j in 1:(i-1)
    @test s[i] != s[j]
    if s[i] == s[j]
        @show (i,j)
    end
end


text:
using Test
using MARY_fRG.Triangulated, MARY_fRG.Basics

@testset "保存相关的功能" begin
    squ1 = Square(Point2D(0., 0.), 2.)
    ltris, ladjs = split_square(squ1, 10)

    save_triangulated(squ1, ltris, ladjs, "squ1")

    @test_broken save_triangulated(ltris[1], ltris, ladjs, "tri1")

    hexa1 = EqHexagon(Point2D(0., 0.), 2.)
    ltris, ladjs = split_hexagon(hexa1, 10)
    save_triangulated(hexa1, ltris, ladjs, "hexa1")
end

@testset "读取相关的功能" begin
    #
    squ1, ltris, ladjs = load_triangulated("squ1")
    save_triangulated(squ1, ltris, ladjs, "squ2")
    orifile = open("squ1.tri", "r")
    cpyfile = open("squ2.tri", "r")
    @test read(orifile) == read(cpyfile)
    rm("squ2.tri")
    orifile = open("squ1.ply", "r")
    cpyfile = open("squ2.ply", "r")
    @test read(orifile) == read(cpyfile)
    rm("squ2.ply")
    orifile = open("squ1.adj", "r")
    cpyfile = open("squ2.adj", "r")
    @test read(orifile) == read(cpyfile)
    rm("squ2.adj")
    #
    hexa1, ltris, ladjs = load_triangulated("hexa1")
    save_triangulated(hexa1, ltris, ladjs, "hexa2")
    orifile = open("hexa1.tri", "r")
    cpyfile = open("hexa2.tri", "r")
    @test read(orifile) == read(cpyfile)
    rm("hexa2.tri")
    orifile = open("hexa1.ply", "r")
    cpyfile = open("hexa2.ply", "r")
    @test read(orifile) == read(cpyfile)
    rm("hexa2.ply")
    orifile = open("hexa1.adj", "r")
    cpyfile = open("hexa2.adj", "r")
    @test read(orifile) == read(cpyfile)
    rm("hexa2.adj")
end

ext: jl


text:
module Crypto

using DiffieHellman: DH
using Pkg.TOML

using ..Keys: Signer
using ..Plugins: Notary, Cypher
using ..Demes: Deme

str(x) = "$x"

sign(data::AbstractString,signer::Signer) = signer.sign(data)
sign(data,signer::Signer) = sign(str(data),signer)

verify(data::AbstractString,signature,deme::Notary) = deme.verify(data,signature)
verify(data,signature,deme::Notary) = verify(str(data),signature,deme)

function hash(x,y,z,notary::Notary)
    inthash = notary.hash("$x $y $z")
    strhash = string(inthash,base=16)
    return Vector{UInt8}(strhash)
end

function wrapsigned(value::BigInt,signer::Signer)
    signature = sign(value,signer)
    signaturedict = Dict(signature)
    dict = Dict("value"=>string(value,base=16),"signature"=>signaturedict)
    io = IOBuffer()
    TOML.print(io,dict)
    return take!(io)
end

function unwrapsigned(envelope::Vector{UInt8},notary::Notary)
    dict = TOML.parse(String(copy(envelope)))
    value = parse(BigInt,dict["value"],base=16)
    signature = notary.Signature(dict["signature"])
    id = verify(value,signature,notary) ### The id is of ID type thus there shall be no problem
    return value, id
end

wrapmsg(value::BigInt) = Vector{UInt8}(string(value,base=16))
unwrapmsg(envelope::Vector{UInt8}) = parse(BigInt,String(copy(envelope)),base=16), nothing

DHsym(cypher::Cypher,notary::Notary,signer::Signer) = DH(value->wrapsigned(value,signer),x->unwrapsigned(x,notary),cypher.G,(x,y,z)->hash(x,y,z,notary),cypher.rng)
DHasym(cypher::Cypher,notary::Notary) = DH(wrapmsg,x->unwrapsigned(x,notary),cypher.G,(x,y,z)->hash(x,y,z,notary),cypher.rng)
DHasym(cypher::Cypher,notary::Notary,signer::Signer) = DH(value->wrapsigned(value,signer),unwrapmsg,cypher.G,(x,y,z)->hash(x,y,z,notary),cypher.rng)

DHsym(deme::Deme,signer::Signer) = DHsym(deme.cypher,deme.notary,signer)
DHasym(deme::Deme) = DHasym(deme.cypher,deme.notary)
DHasym(deme::Deme,signer::Signer) = DHasym(deme.cypher,deme.notary,signer)

export DHsym, DHasym

end


text:
@enum Boundary None=1 Sides=2 Center=3
@enum VertexType Inner=1 Outer=2
@enum InnerBoundaryPos  nnw = 1 nne = 2 nee = 3 see =4 sse = 5 ssw = 6 sww = 7 nww = 8 n = 9 s = 10 w = 11 e = 12
@enum SnappedBoundaryPos nw = 1 n = 2 ne = 3 e = 4 se = 5 s = 6 sw =7 w = 8 c = 9

type BoundaryVertex
  vt::VertexType # needed
  #quadrant::POS
  boundary::InnerBoundaryPos # needed
  #crossing_dir::Bool # positive: crossing in positive direction (left->right, bottom->top)
  vertex::vertex_index #needed
  snapped::Nullable{SnappedBoundaryPos}
end

const BoundaryVertices = Array{Tuple{BoundaryVertex, BoundaryVertex}, 1}

type MeshElement
  triangle_indices::Array{triangle_index, 1}
  boundary_type::Boundary
  boundaries::BoundaryVertices
  center::Nullable{vertex_index}
end

type Triangle
  vertex_indices::Array{vertex_index, 1}
end

# these provide indices into the connectors array for each inner template
# bit patterns give the occupation of the edge mid points (boundaries are
# always occupied)
inner_templates = [
  [[1, 5, 13], [5,9,13]], # 0000
  [[1, 3, 13], [13, 3, 9], [3, 5, 9]], # 0001 - south
  [[1, 11, 13], [1, 11, 5], [5, 11, 9]], # 0010 - north
  [[1, 3, 13], [3, 11, 13], [3, 5, 11], [11, 5, 9]], # 0011 - sn
  [[1, 5, 7], [1, 7, 13], [13, 7, 9]], # 0100 - east
  [[1, 5, 7], [1, 7, 15], [7, 9, 13], [15, 7, 13]], # 0101
  [[1, 5, 7], [7, 9, 11], [1, 11, 13], [1, 7, 11]], # 0110
  [[1, 3, 7], [1, 7, 15], [15, 7, 11], [15, 11, 13], [11, 7, 9]], # 0111
  [[1, 5, 15], [5, 15, 9], [15, 9, 13]], # 1000 - west
  [[1, 3, 15], [3, 5, 9], [15, 9, 13], [3, 9, 15]], # 1001
  [[1, 11, 13], [1, 3, 11], [3, 5, 11], [5, 9, 11]], # 1010
  [[1, 3, 15], [3, 11, 15], [3, 5, 11], [5, 9, 11], [15, 11, 13]], # 1011
  [[1, 3, 13], [3, 5, 7], [7, 9, 13], [3, 7, 13]], # 1100
  [[1, 3, 15], [3, 5, 7], [3, 7, 15], [15, 7, 9], [15, 9, 13]], # 1101
  [[1, 11, 13], [1, 3, 11], [3, 5, 7], [3, 7, 11], [7, 9, 11]], # 1110
  [[1, 3, 15], [3, 5, 7], [15, 3, 11], [15, 11, 13], [3, 7, 11], [7, 9, 11]], # 1111
]

type QuadTreeMesh
  quadtree::QuadTree{MeshElement}
  triangles::Array{Triangle, 1}

  function QuadTreeMesh(bb::GeometryTypes.SimpleRectangle{Float64})
    quadtree = QuadTree{MeshElement}(bb)
    triangles = Array{Triangle, 1}()
    new(quadtree, triangles)
  end

end

"""
    triangulate_leave(qt::QuadTree, elIndex::ElIndex)

Triangulates a leave from `qt`.

# Remarks

* If the leave is not cut by any constraint, it is triangulated according to
  one of six templates
* If the leave contains a constraint vertex, the triangulation is star-shaped
* If the leave is cut but a constraint, it is triangulated according to one of
  18 boundary templates
"""
function triangulate_leave(mesh::QuadTreeMesh, elIndex::ElIndex)

  qt = mesh.quadtree
  qtEl = qt.elements[elIndex]
  @assert(!has_child(qt, elIndex))
  @assert(isnull(qt.values[elIndex]))

  # get vertices of quadtree leave
  # first fill in corner points then half-points
  leave_vertex_indices = zeros(vertex_index, 16)
  leave_vertex_indices[1] = qtEl.bbLeftBottomIndex
  leave_vertex_indices[5] = qtEl.bbRightBottomIndex
  leave_vertex_indices[9] = qtEl.bbRightTopIndex
  leave_vertex_indices[13] = qtEl.bbLeftTopIndex

  shape_index = 0
  neighbour = find_neighbour(qt, elIndex, south)
  if (!isnull(neighbour) && has_child(qt, get(neighbour)))
    shape_index = shape_index + 1
    leave_vertex_indices[3] = get_half_vertex(qt, elIndex, get(neighbour), south)
  end
  neighbour = find_neighbour(qt, elIndex, east)
  if (!isnull(neighbour) && has_child(qt, get(neighbour)))
    shape_index = shape_index + 4
    leave_vertex_indices[7] = get_half_vertex(qt, elIndex, get(neighbour), east)
  end
  neighbour = find_neighbour(qt, elIndex, north)
  if (!isnull(neighbour) && has_child(qt, get(neighbour)))
    shape_index = shape_index + 2
    leave_vertex_indices[11] = get_half_vertex(qt, elIndex, get(neighbour), north)
  end
  neighbour = find_neighbour(qt, elIndex, west)
  if (!isnull(neighbour) && has_child(qt, get(neighbour)))
    shape_index = shape_index + 8
    leave_vertex_indices[15] = get_half_vertex(qt, elIndex, get(neighbour), west)
  end

  # add triangulation templates to triangle list
  templates = inner_templates[shape_index + 1]
  triangle_indices = Array{Int64, 1}()
  for t in templates
    template = [leave_vertex_indices[t[1]], leave_vertex_indices[t[2]], leave_vertex_indices[t[3]]]
    triangle = Triangle(template)
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
  end

  # create new mesh element and add to list
  new_mesh_element = MeshElement(triangle_indices, None, Nullable{vertex_index}(), Nullable{vertex_index}(), Nullable{vertex_index}())
  qt.values[elIndex] = Nullable{QuadTreeMeshes.MeshElement}(new_mesh_element)
end

function get_boundary_vertex_dir(vertex::Point, rt::Point, lb::Point)
  # determine tolerance
  dx, dy = rt - lb

  # shape index bit pattern is wens
  if abs(vertex[2] - rt[2]) < 0.1 * dx
    return 11, 2
  end
  if abs(vertex[2] - lb[2]) < 0.1 * dx
    return 3, 1
  end
  if abs(vertex[1] - lb[1]) < 0.1 * dx
    return 15, 8
  end
  if abs(vertex[1] - rt[1]) < 0.1 * dx
    return 7, 4
  end

  @assert(false)
end

# TODO: corner cases: - constraint along element boundary - vertex at half point - vertex in center

function get_inner_boundary_pos_from_coordinate(vertex::Point, rt::Point, lb::Point)
  # determine tolerance
  dx, dy = rt - lb

  isTop = false
  isLeft = false
  isBottom = false
  isRight = false
  dyt = abs(vertex[2] - rt[2])
  dyb = abs(vertex[2] - lb[2])
  if dyt < 0.1 * dx
    isTop = true
  elseif dyb < 0.1 * dx
    isBottom = true
  end
  dxl = abs(vertex[1] - lb[1])
  dxr = abs(vertex[1] - rt[1])
  if dxl < 0.1 * dx
    isLeft = true
  elseif dxr < 0.1 * dx
    isRight = true
  end
  if isTop
    if dxl < 0.5 * dx
      return nnw
    else
      return nne
    end
  elseif isBottom
    if dxl < 0.5 * dx
      return ssw
    else
      return sse
    end
  elseif isLeft
    if dyt < 0.5 * dy
      return nww
    else
      return sww
    end
  elseif isRight
    if dyt < 0.5 * dy
      return nee
    else
      return see
    end
  else
    @assert(false)
  end
end

function triangulate_boundary_leave(mesh::QuadTreeMesh, elIndex::ElIndex, bnd1Index::vertex_index, bnd2Index::vertex_index)
  qt = mesh.quadtree
  qtEl = qt.elements[elIndex]
  @assert(!has_child(qt, elIndex))
  @assert(isnull(qt.values[elIndex]))

  # get vertices of quadtree leave
  # fill in corner points
  leave_vertex_indices = zeros(vertex_index, 18)
  leave_vertex_indices[1] = qtEl.bbLeftBottomIndex
  leave_vertex_indices[5] = qtEl.bbRightBottomIndex
  leave_vertex_indices[9] = qtEl.bbRightTopIndex
  leave_vertex_indices[13] = qtEl.bbLeftTopIndex
  leave_vertex_indices[17] = bnd1Index
  leave_vertex_indices[18] = bnd2Index

  # boundary leaves cannot have half-grid points
  neighbour = find_neighbour(qt, elIndex, south)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))
  neighbour = find_neighbour(qt, elIndex, north)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))
  neighbour = find_neighbour(qt, elIndex, west)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))
  neighbour = find_neighbour(qt, elIndex, east)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))

  rt = qt.vertices[qtEl.bbRightTopIndex]
  lb = qt.vertices[qtEl.bbLeftBottomIndex]

  # we pretend bnd1 indices are at half grid points -
  # this only concerns the topology
  bnd1Pos, bnd1Shape = get_boundary_vertex_dir(qt.vertices[bnd1Index], rt, lb)
  bnd2Pos, bnd2Shape = get_boundary_vertex_dir(qt.vertices[bnd2Index], rt, lb)
  @assert(bnd1Shape != bnd2Shape) # TODO: this is a special corner case..

  leave_vertex_indices[bnd1Pos] = bnd1Index
  leave_vertex_indices[bnd2Pos] = bnd2Index
  shape_index = bnd1Shape + bnd2Shape

  # add triangulation templates to triangle list
  templates = inner_templates[shape_index + 1]
  triangle_indices = Array{Int64, 1}()
  for t in templates
    template = [leave_vertex_indices[t[1]], leave_vertex_indices[t[2]], leave_vertex_indices[t[3]]]
    triangle = Triangle(template)
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
  end

  # we need to classify the boundary vertices from that triangulation
  # we do not check if the given vertex is actually *on* the boundary
  # we just assume it is
  pos1 = get_inner_boundary_pos_from_coordinate(qt.vertices[bnd1Index], rt, lb)
  pos2 = get_inner_boundary_pos_from_coordinate(qt.vertices[bnd2Index], rt, lb)
  bv1 = BoundaryVertex(Outer, pos1, bnd1Index, Nullable{SnappedBoundaryVertex}())
  bv2 = BoundaryVertex(Outer, pos2, bnd2Index, Nullable{SnappedBoundaryVertex}())

  # create new mesh element and add to list
  new_mesh_element = MeshElement(triangle_indices, Sides, [(bv1, bv2)], Nullable{vertex_index}())
  qt.values[elIndex] = Nullable{MeshElement}(new_mesh_element)
end

function get_vertex_pos_on_boundary(cp::Point, bd::Point)
  if bd[1] < cp[1]
    if bd[2] < cp[2]
      return northWest
    else
      return southWest
    end
  else
    if bd[2] < cp[2]
      return northEast
    else
      return southEast
    end
  end
end

function triangulate_boundary_leave_with_vertex(mesh::QuadTreeMesh, elIndex::ElIndex, bnd1Index::vertex_index, bnd2Index::vertex_index, vIndex::vertex_index)
  qt = mesh.quadtree
  qtEl = qt.elements[elIndex]
  @assert(!has_child(qt, elIndex))
  @assert(isnull(qt.values[elIndex]))

  # TODO:
  # - get center point and determine position of boundary points
  # - create triangles nw - v - ne or nw - v - b/v - ne - b
  #   and accordingly for all other segments
  # TODO: what about two boundary vertices on same element? disallow for now.

  # TODO: remove for production
  # boundary leaves cannot have half-grid points
  neighbour = find_neighbour(qt, elIndex, south)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))
  neighbour = find_neighbour(qt, elIndex, north)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))
  neighbour = find_neighbour(qt, elIndex, west)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))
  neighbour = find_neighbour(qt, elIndex, east)
  assert(isnull(neighbour) || !has_child(qt, get(neighbour)))

  # get center point and determine position of boundary leave
  rt = qt.vertices[qtEl.bbRightTopIndex]
  lb = qt.vertices[qtEl.bbLeftBottomIndex]
  c = 0.5 * (rt + lb)
  bd1pos = get_vertex_pos_on_boundary(c, qt.vertices[bnd1Index])
  bd2pos = get_vertex_pos_on_boundary(c, qt.vertices[bnd2Index])

  triangle_indices = Array{Int64, 1}()
  if bd1pos != northWest && bd2pos != northWest
      triangle = Triangle([qtEl.bbLeftTopIndex, vIndex, qtEl.bbRightTopIndex])
      push!(mesh.triangles, triangle)
      push!(triangle_indices, length(mesh.triangles))
  else
    if bd1pos == northWest
      bdIndex = bnd1Index
    else
      bdIndex = bnd2Index
    end
    triangle = Triangle([qtEl.bbLeftTopIndex, vIndex, bdIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
    triangle = Triangle([qtEl.bbRightTopIndex, bdIndex, vIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
  end
  if bd1pos != northEast && bd2pos != northEast
      triangle = Triangle([qtEl.bbLeftTopIndex, vIndex, qtEl.bbRightTopIndex])
      push!(mesh.triangles, triangle)
      push!(triangle_indices, length(mesh.triangles))
  else
    if bd1pos == northEast
      bdIndex = bnd1Index
    else
      bdIndex = bnd2Index
    end
    triangle = Triangle([qtEl.bbLeftTopIndex, vIndex, bdIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
    triangle = Triangle([qtEl.bbRightTopIndex, bdIndex, vIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
  end
  if bd1pos != southWest && bd2pos != southWest
      triangle = Triangle([qtEl.bbLeftBottomIndex, vIndex, qtEl.bbRightBottomIndex])
      push!(mesh.triangles, triangle)
      push!(triangle_indices, length(mesh.triangles))
  else
    if bd1pos == southWest
      bdIndex = bnd1Index
    else
      bdIndex = bnd2Index
    end
    triangle = Triangle([qtEl.bbLeftBottomIndex, vIndex, bdIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
    triangle = Triangle([qtEl.bbRightBottomIndex, bdIndex, vIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
  end
  if bd1pos != southEast && bd2pos != southEast
      triangle = Triangle([qtEl.bbLeftBottomIndex, vIndex, qtEl.bbRightBottomIndex])
      push!(mesh.triangles, triangle)
      push!(triangle_indices, length(mesh.triangles))
  else
    if bd1pos == southEast
      bdIndex = bnd1Index
    else
      bdIndex = bnd2Index
    end
    triangle = Triangle([qtEl.bbLeftBottomIndex, vIndex, bdIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
    triangle = Triangle([qtEl.bbRightBottomIndex, bdIndex, vIndex])
    push!(mesh.triangles, triangle)
    push!(triangle_indices, length(mesh.triangles))
  end

  # create new mesh element and add to list
  nbvv = BoundaryVertex(Inner, nnw, vIndex, Nullable{SnappedBoundaryVertex}())
  pos1 = get_inner_boundary_pos_from_coordinate(qt.vertices[bnd1Index], rt, lb)
  nbv1 = BoundaryVertex(Outer, pos1, bnd1Index, Nullable{SnappedBoundaryVertex}())
  pos2 = get_inner_boundary_pos_from_coordinate(qt.vertices[bnd2Index], rt, lb)
  nbv2 = BoundaryVertex(Outer, pos2, bnd2Index, Nullable{SnappedBoundaryVertex}())

  new_mesh_element = MeshElement(triangle_indices, Center, [(nbv1, nbvv), (nbvv, nbv2)], Nullable{vertex_index}(vIndex))
  qt.values[elIndex] = Nullable{QuadTreeMeshes.MeshElement}(new_mesh_element)
end

type LeaveBoundarySegments
  nnw_segment::GeometryTypes.Simplex{2, Point}
  nne_segment::GeometryTypes.Simplex{2, Point}
  nee_segment::GeometryTypes.Simplex{2, Point}
  see_segment::GeometryTypes.Simplex{2, Point}
  sse_segment::GeometryTypes.Simplex{2, Point}
  ssw_segment::GeometryTypes.Simplex{2, Point}
  sww_segment::GeometryTypes.Simplex{2, Point}
  nww_segment::GeometryTypes.Simplex{2, Point}
  n_segment::GeometryTypes.Simplex{2, Point}
  s_segment::GeometryTypes.Simplex{2, Point}
  w_segment::GeometryTypes.Simplex{2, Point}
  e_segment::GeometryTypes.Simplex{2, Point}

  function LeaveBoundarySegments(qt::QuadTree, elIndex::ElIndex)
    qt_element = qt.elements[elIndex]

    # new children
    nw_element = qt.elements[get(qt_element.northWest)]
    ne_element = qt.elements[get(qt_element.northEast)]
    sw_element = qt.elements[get(qt_element.southWest)]
    se_element = qt.elements[get(qt_element.southEast)]

    # all relevant vertices
    nw_vertex = qt_element.bbLeftTopIndex
    n_vertex = nw_element.bbRightTopIndex
    ne_vertex = ne_element.bbRightTopIndex
    e_vertex = ne_element.bbRightBottomIndex
    se_vertex = se_element.bbRightBottomIndex
    s_vertex = se_element.bbLeftBottomIndex
    sw_vertex = sw_element.bbLeftBottomIndex
    w_vertex = sw_element.bbLeftTopIndex
    c_vertex = nw_element.bbRightBottomIndex

    # create segments - make sure they always point in positive x/y dir
    nnw_segment = GeometryTypes.LineSegment(qt.vertices[nw_vertex], qt.vertices[n_vertex])
    nne_segment = GeometryTypes.LineSegment(qt.vertices[n_vertex], qt.vertices[ne_vertex])
    nee_segment = GeometryTypes.LineSegment(qt.vertices[e_vertex], qt.vertices[ne_vertex])
    see_segment = GeometryTypes.LineSegment(qt.vertices[se_vertex], qt.vertices[e_vertex])
    sse_segment = GeometryTypes.LineSegment(qt.vertices[s_vertex], qt.vertices[se_vertex])
    ssw_segment = GeometryTypes.LineSegment(qt.vertices[sw_vertex], qt.vertices[s_vertex])
    sww_segment = GeometryTypes.LineSegment(qt.vertices[sw_vertex], qt.vertices[w_vertex])
    nww_segment = GeometryTypes.LineSegment(qt.vertices[w_vertex], qt.vertices[nw_vertex])
    n_segment = GeometryTypes.LineSegment(qt.vertices[c_vertex], qt.vertices[n_vertex])
    s_segment = GeometryTypes.LineSegment(qt.vertices[s_vertex], qt.vertices[c_vertex])
    w_segment = GeometryTypes.LineSegment(qt.vertices[w_vertex], qt.vertices[c_vertex])
    e_segment = GeometryTypes.LineSegment(qt.vertices[c_vertex], qt.vertices[e_vertex])

    new(nnw_segment, nne_segment, nee_segment, see_segment, sse_segment, ssw_segment, sww_segment, nww_segment, n_segment, s_segment, w_segment, e_segment)
  end
end

function String(x::LeaveBoundarySegments)
  restr = """
$(x.nnw_segment)
$(x.nne_segment)
$(x.nee_segment)
$(x.see_segment)
$(x.sse_segment)
$(x.ssw_segment)
$(x.sww_segment)
$(x.nww_segment)
$(x.nww_segment)
$(x.n_segment)
$(x.e_segment)
$(x.s_segment)
$(x.w_segment)
"""
  return restr
end

function intersect(ls1::GeometryTypes.LineSegment, ls2::GeometryTypes.LineSegment)
  # parametric solution - https://gamedev.stackexchange.com/questions/44720/line-intersection-from-parametric-equation
  # returns parameter of first line segment
  s1, e1 = ls1
  s2, e2 = ls2
  s1x, s1y = s1
  s2x, s2y = s2
  e1x, e1y = e1
  e2x, e2y = e2
  # ls1 = a + t b (0 <= t <= 1)
  ax=s1x
  ay=s1y
  b = e1 - s1
  bx, by = b
  cx=s2x
  cy=s2y
  d = e2 - s2
  dx, dy = d
  if ((dx*by-dy*bx) == 0) || ((bx*dy-by*dx) == 0)
    return false, 0, Point(0.,0.), false
  end

  u=(bx*(cy-ay) + by*(ax-cx))/(dx*by-dy*bx)
  t=(dx*(ay-cy) + dy*(cx-ax))/(bx*dy-by*dx)

  if u >= 0 && u <= 1 && t >= 0 && t <= 1
    # get intersection point and intersection dir
    p = Point(ax + t * bx, ay + t * by)
    s = (dx != 0 && bx > 0) || (dy != 0 && by > 0)
    #TODO: I don't think we need the direction
    return true, t, p, s
  else
    return false, 0, Point(0.,0.), true
  end
end

function check_for_single_endpoint(ls::GeometryTypes.LineSegment, isp::Point, bpl::SnappedBoundaryPos, ibpr::SnappedBoundaryPos)
  s, e = ls

  ds = dot((s-isp), (s-isp))
  if (ds < 1e-5^2)
    return Nullable{SnappedBoundaryPos}(ibpl)
  end
  de = dot((e-isp), (e-isp))
  if (de < 1e-5^2)
    return Nullable{SnappedBoundaryPos}(ibpr)
  end

  return Nullable{SnappedBoundaryPos}()
end

function fill_if_intersects(ls1::GeometryTypes.LineSegment, ls2::GeometryTypes.LineSegment, isp::Array{Tuple{Float64, Point, Bool, InnerBoundaryPos, SnappedBoundaryPos}, 1}, ibp::InnerBoundaryPos, sp1::SnappedBoundaryPos, sp2::SnappedBoundaryPos)
  i, t, p, s = intersect(ls1, ls2)
  if i
    #print("Pushing point $p from $ls1 -> $ls2\n")
    # snap boundary points onto endpoints if they are close enough
    snapped = check_for_single_endpoint(ls2, i, sp1, sp2)
    push!(isp, (t, p, s, ibp, snapped))
  end
end

function intersect_with_leave(lbs::LeaveBoundarySegments, ls::GeometryTypes.Simplex{2, Point})
  # get all intersections with all relevant boundary segments
  isp = Array{Tuple{Float64, Point, Bool, InnerBoundaryPos, SnappedBoundaryPos}, 1}()
  fill_if_intersects(ls, lbs.nnw_segment, isp, nnw, nw, n)
  fill_if_intersects(ls, lbs.nne_segment, isp, nne, n, ne)
  fill_if_intersects(ls, lbs.nee_segment, isp, nee, e, ne)
  fill_if_intersects(ls, lbs.see_segment, isp, see, se, e)
  fill_if_intersects(ls, lbs.sse_segment, isp, sse, s, se)
  fill_if_intersects(ls, lbs.ssw_segment, isp, ssw, sw, s)
  fill_if_intersects(ls, lbs.sww_segment, isp, sww, sw, w)
  fill_if_intersects(ls, lbs.nww_segment, isp, nww, w, nw)
  fill_if_intersects(ls, lbs.n_segment, isp, n, c, n)
  fill_if_intersects(ls, lbs.s_segment, isp, s, s, c)
  fill_if_intersects(ls, lbs.w_segment, isp, w, w, c)
  fill_if_intersects(ls, lbs.e_segment, isp, e, c, e)

  # sort them according to distance from start point
  sort!(isp, lt = (is1, is2) -> return is1[1] < is2[1])

  return isp
end

function get_quadrant_from_boundary_position(qt::QuadTree, elIndex::ElIndex, intersectionPoint::BoundaryVertex)
  if intersectionPoint.vt == Outer
    ips = intersectionPoint.boundary
    print_with_color(:blue, "Checking intersection point $ips.\n")
    retVal = Array{POS, 1}()
    if !isnull(intersectionPoint.snapped)
      snapped = get(intersectionPoint.snapped)
      if snapped == nw || snapped == n || snapped == w || snapped == c
        push!(retVal, northWest)
      end
      if snapped == n || snapped == ne || snapped == e || snapped == c
        push!(retVal, northEast)
      end
      if snapped == c || snapped == e || snapped == s || snapped == se
        push!(retVal, southEast)
      end
      if snapped == c || snapped == w || snapped == s || snapped == sw
        push!(retVal, southWest)
      end
    else
      if ips == nnw || ips == nww || ips == n || ips == w
        push!(retVal, northWest)
      end
      if ips == nne || ips == nee || ips == n || ips == e
        push!(retVal, northEast)
      end
      if ips == sse || ips == see || ips == s || ips == e
        push!(retVal, southEast)
      end
      if ips == ssw || ips == sww || ips == s || ips == w
        push!(retVal, southWest)
      end
    end
    return retVal
  elseif intersectionPoint.vt == Inner
    vpp = qt.vertices[intersectionPoint.vertex]
    pos = get_child_positions_from_coordinate(qt, elIndex, vpp)
    return pos
  end
end

function find_common_positions(first::Array{POS, 1}, second::Array{POS,1})
  # we assume the input is sorted (which is the case if it comes from get_quadrant_from_boundary_position)
  firstIdx = 1
  secondIdx = 1
  common_positions = Array{POS, 1}
  while firstIdx <= length(first) && secondIdx <= length(second)
    if first[firstIdx] == second[secondIdx]
      push!(common_positions, first[firstIdx])
    elseif first[firstIdx] < second[secondIdx] # apparently we can compare enums..
      firstIdx = firstIdx + 1
    else
      secondIdx = secondIdx + 1
    end
  end
  @assert(length(common_positions)>=1, "No common position found for $first and $second.")
  return common_positions
end

# we assume that the members of isp have corrrect boundary, vertex_type, crossing_dir
function forward_boundaries_to_leaves(qt::QuadTree, parent_element::ElIndex, isp::Array{BoundaryVertex})
    #      * for first intersection: assert that second intersection is in same quadrant and fill in structure in child
    #      * for second intersection: assert that third intersection is in same quadrant and fill in structure in child
    #      * etc.. until last pair was found
    @assert(length(isp) >= 2)
    curSecondIdx = 2
    firstIsp = isp[1]
    while curSecondIdx <= length(isp)
      secondIsp = isp[curSecondIdx]

      if firstIsp.vertex != secondIsp.vertex
        # intersection points can belong to up to four quadrants (if it's on an inner boundary)
        # find all common positions
        print_with_color(:green, "Testing  $(qt.vertices[firstIsp.vertex]) and $(qt.vertices[secondIsp.vertex]).")
        firstPoss = get_quadrant_from_boundary_position(qt, parent_element, firstIsp)
        secondPoss = get_quadrant_from_boundary_position(qt, parent_element, secondIsp)
        poss = find_common_positions(firstPoss, secondPoss)

        # for each common quadrant create boundary vertices from intersection points
        for pos in poss
          leaveIndex = get_leave_from_pos(qt, parent_element, pos)
          leave = get(qt.values[leaveIndex])
          bv1 = BoundaryVertex(firstIsp.vt, firstIsp.boundary, firstIsp.vertex, Nullable{SnappedBoundaryVertex}())
          bv2 = BoundaryVertex(secondIsp.vt, secondIsp.boundary, secondIsp.vertex, Nullable{SnappedBoundaryVertex}())
          push!(leave.boundaries, (bv1, bv2))
          if (firstIsp.vt == Inner)
            leave.center = Nullable{vertex_index}(firstIsp.vertex)
            leave.boundary_type = Center
          elseif  secondIsp.vt == Inner
            leave.boundary_type = Center
            leave.center = Nullable{vertex_index}(secondIsp.vertex)
          else
            leave.boundary_type = Sides
          end
        end
      end

      firstIsp = secondIsp
      curSecondIdx += 1
    end
end

function get_snapped_vertex(qt::QuadTree, parent_element::ElIndex, snapped::SnappedBoundaryVertex)
  qtEl = qt.elements[parent_element]
  if snapped == nw
    return qtEl.bbLeftTopIndex
  elseif snapped == ne
    return qtEl.bbRightTopIndex
  elseif snapped == sw
    return qtEl.bbLeftBottomIndex
  elseif snapped == se
    return qtEl.bbRightBottomIndex
  elseif snapped == n
    northWestEl = qt.elements[get(qtEl.northWest)]
    return northWestEl.bbRightTopIndex
  elseif snapped == e
    northEastEl == qt.elements[get(qtEl.northEast)]
    return northEastEl.bbRightBottomIndex
  elseif snapped == s
    southEastEl == qt.elements[get(qtEl.southEast)]
    return southEastEl.bbLeftBottomIndex
  elseif snapped == w
    southWestEl == qt.elements[get(qtEl.southWest)]
    return southWestEl.bbLeftTopIndex
  elseif snapped == c
    southWestEl == qt.elements[get(qtEl.southWest)]
    return southWestEl.bbRightTopIndex
  end
end

function push_or_get_snapped_vertex(qt::QuadTree, parent_element::ElIndex, p::Point, snapped::Nullable{SnappedBoundaryVertex})
  if isnull(snapped)
    push!(qt.vertices, p)
    return length(qt.vertices)
  else
    return get_snapped_vertex(qt, parent_element, get(snapped))
  end
end

function propagate_intersections(qt::QuadTree, parent_element::ElIndex, bndy::Tuple{BoundaryVertex, BoundaryVertex})
  # comppute intersection points of boundary with all segments of subdivided leave
  # TODO: can probably pull that out of the function
  lbs = LeaveBoundarySegments(qt, parent_element)
  b1, b2 = bndy
  ls1p1 = qt.vertices[b1.vertex]
  ls1p2 = qt.vertices[b2.vertex]
  ls1 = GeometryTypes.Simplex{2, Point}(ls1p1, ls1p2)
  isps = intersect_with_leave(lbs, ls1)

  # transform intersection points into BoundaryVertex
  allBoundaries = Array{BoundaryVertex, 1}()
  if b1.vt == Inner
    startIndex = 1
    push!(allBoundaries, b1)
  else
    # ignore first intersection point - it is a bndy point of the parent
    startIndex = 2
    push!(allBoundaries, b1)
  end

  # create vertices for all intersections except
  # first one (if it's an inner boundary) and last
  # one
  usedIsps = isps[startIndex:length(isps)-1]
  for isp in usedIsps
    # create vertex and make it into a boundary vertex
    newVertexIndex = push_or_get_snapped_vertex(qt, parent_element, isp[2], isp[5])
    nbv = BoundaryVertex(Outer, isp[4], newVertexIndex, Nullable{SnappedBoundaryVertex}())
    push!(allBoundaries,nbv)
  end

  # if second boundary is an inner vertex, add last intersection as vertex and push inner vertex
  # otherwise just push second vertex
  if b2.vt == Inner
    # push last intersection and end with inner
    if length(isps) >= startIndex
      isp = isps[length(isps)]
      newVertexIndex = push_or_get_snapped_vertex(qt, parent_element, isp[2], isp[5])
      nbv = BoundaryVertex(Outer, isp[4], newVertexIndex, Nullable{SnappedBoundaryVertex}())
      push!(allBoundaries,nbv)
    end
    push!(allBoundaries, b2)
  else
    push!(allBoundaries, b2)
  end

  for b in allBoundaries
    print_with_color(:green, "$(b.vertex) => $(qt.vertices[b.vertex])\n")
  end
  for i in isps
    print_with_color(:blue, "$(i[2])\n")
  end

  forward_boundaries_to_leaves(qt, parent_element, allBoundaries)
end

"""
  OnChildrenCreated(qt::QuadTree, elIndex::ElIndex)

Will be called when quadtree element `elIndex` is subdivided. This method examines
the element for boundaries and constraint vertices and propagates them to the
children if necessary.

# Remarks
* A constraint vertex is moved into the corresponding child element
* If a boundary is present, it is intersected with the bounding boxes of the inner
  elements and the intersection vertices are stored in the child elements.
"""
function OnChildrenCreated(qt::QuadTree, elIndex::ElIndex)

  mesh_element = get(qt.values[elIndex])
  qt_element = qt.elements[elIndex]

  # create empty leaves for all children
  qt.values[get(qt_element.northWest)] = MeshElement(Array{triangle_index, 1}(), None, BoundaryVertices(), Nullable{vertex_index}())
  qt.values[get(qt_element.northEast)] = MeshElement(Array{triangle_index, 1}(), None, BoundaryVertices(), Nullable{vertex_index}())
  qt.values[get(qt_element.southWest)] = MeshElement(Array{triangle_index, 1}(), None, BoundaryVertices(), Nullable{vertex_index}())
  qt.values[get(qt_element.southEast)] = MeshElement(Array{triangle_index, 1}(), None, BoundaryVertices(), Nullable{vertex_index}())

  for bnd in mesh_element.boundaries
    propagate_intersections(qt, elIndex, bnd)
  end
end


text:
using SimpleTest
using Test

@testset "SimpleTest.jl" begin
    @test simple_operation(2.0, 3.0) == 5
    @test simple_operation(2.1, 3.1)==5.2
    @test simple_operation(21.0, 31.0)== 52.0
end


text:
using PkgTemplates
using Test
using Dates
using LibGit2
using Pkg

import PkgTemplates: badges, version_floor, substitute, read_license, gen_file, gen_readme,
    gen_tests, gen_license, gen_gitignore, gen_plugin, show_license, LICENSES,
    LICENSE_DIR, Plugin, GenericPlugin, CustomPlugin, Badge, format, interactive,
    DEFAULTS_DIR, Documenter

mktempdir() do temp_dir
    mkdir(joinpath(temp_dir, "dev"))
    pushfirst!(DEPOT_PATH, temp_dir)
    cd(temp_dir) do
        @testset "PkgTemplates.jl" begin
            include("tests.jl")
        end
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the API to convert from small Euler angles to rotation
#   representations.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export smallangle_to_rot

"""
    smallangle_to_rot([T,] θx::Number, θy::Number, θz::Number[; normalize = true])

Create a rotation description of type `T` from three small rotations of angles
`θx`, `θy`, and `θz` [rad] about the axes X, Y, and Z, respectively.

The type `T` of the rotation description can be `DCM` or `Quaternion`. If the
type `T` is not specified, then if defaults to `DCM`.

If `T` is `DCM`, then the resulting matrix will be orthonormalized using the
`orthonormalize` function if the keyword `normalize` is `true`.

# Example

```jldoctest
julia> dcm = smallangle_to_rot(+0.01, -0.01, +0.01)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.9999     0.00989903  0.010098
 -0.009999   0.999901    0.00989802
 -0.009999  -0.009998    0.9999

julia> dcm = smallangle_to_rot(+0.01, -0.01, +0.01; normalize = false)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  1.0    0.01  0.01
 -0.01   1.0   0.01
 -0.01  -0.01  1.0

julia> q = smallangle_to_rot(Quaternion, +0.01, -0.01, +0.01)
Quaternion{Float64}:
  + 0.999963 + 0.00499981⋅i - 0.00499981⋅j + 0.00499981⋅k
```
"""
@inline function smallangle_to_rot(
    θx::Number,
    θy::Number,
    θz::Number;
    normalize = true
)
    return smallangle_to_dcm(θx, θy, θz; normalize = normalize)
end

@inline function smallangle_to_rot(
    ::Type{DCM},
    θx::Number,
    θy::Number,
    θz::Number;
    normalize = true
)
    return smallangle_to_dcm(θx, θy, θz; normalize = normalize)
end

@inline function smallangle_to_rot(
    ::Type{Quaternion},
    θx::Number,
    θy::Number,
    θz::Number
)
    return smallangle_to_quat(θx, θy, θz)
end


text:
using Revise
using FastGroupBy, BenchmarkTools

#const N = Int64(2e9/8)
#const N = 100_000_000
#const N = UInt32(2^31)
const N = UInt32(2^30)

#const N = Int(2^31-1)
const K = UInt32(100)
#const id4 = rand(1:K, N)
const id6 = rand(Int32(1):Int32(round(N/K)), N)
const v1 =  rand(Int32(1):Int32(5), N)
#const v1 = similar(id6)
#const v3 =  rand(round.(rand(100)*100,4), N)

#@time sumby(id6, v1)

# 526 seconds for 2billion
#using SortingAlgorithms
#@time sumby(id6, v1)

# @belapsed sumby(id4,v1) # 4.8
# @belapsed sumby(id4,v3) #4.7
#@time psumby(id4,v1); # 7.8
#@time psumby(id4,v3); #7.7

#
# @belapsed sumby(id6,v1) #64.5
# @belapsed sumby(id6,v3) #68.3
#@belapsed psumby(id6,v1) #28.3
#@belapsed psumby(id6,v3) #37.2


text:
#
# Testing routines: computing the lists with naive algorithms
#

# For a single set of molecules
function naive_md(x, n_atoms_per_molecule_x::Int, box::Box)
    x_list = init_list(x, i -> mol_index(i, n_atoms_per_molecule_x))
    for i in 1:length(x)-1
        vx = x[i]
        for j in i+1:length(x)
            vy = x[j]
            imol = mol_index(i, n_atoms_per_molecule_x)
            jmol = mol_index(j, n_atoms_per_molecule_x)
            if imol == jmol
                continue
            end
            vy_wrapped = CellListMap.wrap_relative_to(vy, vx, box)
            d = sqrt(sum(abs2, vy_wrapped - vx))
            if d > box.cutoff
                continue
            end
            if d < x_list[imol].d
                x_list[imol] = MinimumDistance(true, i, j, d)
            end
            if d < x_list[jmol].d
                x_list[jmol] = MinimumDistance(true, j, i, d)
            end
        end
    end
    return x_list
end

# For disjoint sets, returning only one list
function naive_md(x, y, n_atoms_per_molecule_x::Int, box::Box)
    x_list = init_list(x, i -> mol_index(i, n_atoms_per_molecule_x))
    for (i, vx) in pairs(x)
        for (j, vy) in pairs(y)
            vy_wrapped = CellListMap.wrap_relative_to(vy, vx, box)
            d = sqrt(sum(abs2, vy_wrapped - vx))
            if d > box.cutoff
                continue
            end
            imol = mol_index(i, n_atoms_per_molecule_x)
            if d < x_list[imol].d
                x_list[imol] = MinimumDistance(true, i, j, d)
            end
        end
    end
    return x_list
end

# For disjoint sets, returning both lists
function naive_md(x, y, n_atoms_per_molecule_x::Int, n_atoms_per_molecule_y::Int, box::Box)
    x_list = init_list(x, i -> mol_index(i, n_atoms_per_molecule_x))
    y_list = init_list(y, i -> mol_index(i, n_atoms_per_molecule_y))
    for (i, vx) in pairs(x)
        for (j, vy) in pairs(y)
            vy_wrapped = CellListMap.wrap_relative_to(vy, vx, box)
            d = sqrt(sum(abs2, vy_wrapped - vx))
            if d > box.cutoff
                continue
            end
            imol = mol_index(i, n_atoms_per_molecule_x)
            if d < x_list[imol].d
                x_list[imol] = MinimumDistance(true, i, j, d)
            end
            jmol = mol_index(j, n_atoms_per_molecule_y)
            if d < y_list[jmol].d
                y_list[jmol] = MinimumDistance(true, j, i, d)
            end
        end
    end
    return x_list, y_list
end

#
# Comparison functions
#
import Base.isapprox
function isapprox(
    list1::List,
    list2::List,
)
    length(list1) != length(list2) && return false
    for i in eachindex(list1)
        list1[i].within_cutoff != list2[i].within_cutoff && return false
        list1[i].i != list2[i].i && return false
        list1[i].j != list2[i].j && return false
        !(list1[i].d ≈ list2[i].d) && return false
    end
    return true
end
function isapprox(
    lists1::Tuple{T,T},
    lists2::Tuple{T,T},
) where {T<:List}
    !(lists1[1] ≈ lists2[1]) && return false
    !(lists1[2] ≈ lists2[2]) && return false
    return true
end

#
# Auxiliary functions to plot a 2D molecular arrangement for
# inspection and visualization
#
function plot_mol!(
    p,
    x;
    cycle=false,
    markercolor=nothing
)
    for i in 1:length(x)-1
        xat = [ x[i][1], x[i+1][1] ]
        yat = [ x[i][2], x[i+1][2] ]
        Main.plot!(p, xat, yat, label=:none, linewidth=2, color=:black)
    end
    if cycle
        xat = [ x[end][1], x[1][1] ]
        yat = [ x[end][2], x[1][2] ]
        Main.plot!(p, xat, yat, label=:none, linewidth=2, color=:black)
    end
    if isnothing(markercolor)
        Main.scatter!(p, Tuple.(x), label=:none, markersize=4.0)
    else
        Main.scatter!(p, Tuple.(x), label=:none, markersize=4.0, markercolor=markercolor)
    end
    return p
end

function plot_mol!(
    p,
    x,
    n_atoms_per_molecule::Int;
    cycle=false,
    markercolor=nothing
)
    for i in 1:n_atoms_per_molecule:length(x)
        plot_mol!(
            p,
            @view(x[i:i+n_atoms_per_molecule-1]);
            cycle=cycle,
            markercolor=markercolor
        )
        if isnothing(markercolor)
            markercolor = p.series_list[end][:markercolor]
        end
    end
    return p
end

function plot_md!(
    p,
    x,
    n_atoms_per_molecule_x::Int,
    y,
    n_atoms_per_molecule_y::Int,
    md::List;
    x_cycle=false,
    y_cycle=false
)
    for pair in md
        !pair.within_cutoff && continue
        xp = [ x[pair.i][1], y[pair.j][1] ]
        yp = [ x[pair.i][2], y[pair.j][2] ]
        Main.plot!(p, xp, yp, label=:none, linewidth=2, color=:black, alpha=0.3)
    end
    plot_mol!(p, x, n_atoms_per_molecule_x; cycle=x_cycle, markercolor=:red)
    plot_mol!(p, y, n_atoms_per_molecule_y; cycle=y_cycle, markercolor=:blue)
end

download_example() =
    Main.PDBTools.readPDB(download("https://raw.githubusercontent.com/m3g/ComplexMixtures.jl/master/test/data/NAMD/structure.pdb"))

inds(i) = mol_index(i,3)
function iterate_lists(nsteps, list, list_threaded, water, protein, box, cl, aux; parallel = true)
   for i in 1:nsteps
       # water and protein coordinates, and box, could change here
       cl = UpdateCellList!(water, protein, box, cl, aux; parallel = parallel)
       minimum_distances!(inds, list, box, cl;
           parallel = parallel,
           list_threaded = list_threaded
       )
       # list was updated
   end
end


text:
function hubbard_1d(; N::Int, t=1.0, U=0.0)
  ampo = OpSum()
  for b in 1:(N - 1)
    ampo .+= -t, "Cdagup", b, "Cup", b + 1
    ampo .+= -t, "Cdagup", b + 1, "Cup", b
    ampo .+= -t, "Cdagdn", b, "Cdn", b + 1
    ampo .+= -t, "Cdagdn", b + 1, "Cdn", b
  end
  if U ≠ 0
    for n in 1:N
      ampo .+= U, "Nupdn", n
    end
  end
  return ampo
end

function hubbard_2d(; Nx::Int, Ny::Int, t=1.0, U=0.0, yperiodic::Bool=true)
  N = Nx * Ny
  lattice = square_lattice(Nx, Ny; yperiodic=yperiodic)
  ampo = OpSum()
  for b in lattice
    ampo .+= -t, "Cdagup", b.s1, "Cup", b.s2
    ampo .+= -t, "Cdagup", b.s2, "Cup", b.s1
    ampo .+= -t, "Cdagdn", b.s1, "Cdn", b.s2
    ampo .+= -t, "Cdagdn", b.s2, "Cdn", b.s1
  end
  if U ≠ 0
    for n in 1:N
      ampo .+= U, "Nupdn", n
    end
  end
  return ampo
end

function hubbard_2d_ky(; Nx::Int, Ny::Int, t=1.0, U=0.0)
  ampo = OpSum()
  for x in 0:(Nx - 1)
    for ky in 0:(Ny - 1)
      s = x * Ny + ky + 1
      disp = -2 * t * cos((2 * π / Ny) * ky)
      if abs(disp) > 1e-12
        ampo .+= disp, "Nup", s
        ampo .+= disp, "Ndn", s
      end
    end
  end
  for x in 0:(Nx - 2)
    for ky in 0:(Ny - 1)
      s1 = x * Ny + ky + 1
      s2 = (x + 1) * Ny + ky + 1
      ampo .+= -t, "Cdagup", s1, "Cup", s2
      ampo .+= -t, "Cdagup", s2, "Cup", s1
      ampo .+= -t, "Cdagdn", s1, "Cdn", s2
      ampo .+= -t, "Cdagdn", s2, "Cdn", s1
    end
  end
  if U ≠ 0
    for x in 0:(Nx - 1)
      for ky in 0:(Ny - 1)
        for py in 0:(Ny - 1)
          for qy in 0:(Ny - 1)
            s1 = x * Ny + (ky + qy + Ny) % Ny + 1
            s2 = x * Ny + (py - qy + Ny) % Ny + 1
            s3 = x * Ny + py + 1
            s4 = x * Ny + ky + 1
            ampo .+= (U / Ny), "Cdagdn", s1, "Cdagup", s2, "Cup", s3, "Cdn", s4
          end
        end
      end
    end
  end
  return ampo
end

function hubbard(; Nx::Int, Ny::Int=1, t=1.0, U=0.0, yperiodic::Bool=true, ky::Bool=false)
  if Ny == 1
    ampo = hubbard_1d(; N=Nx, t=t, U=U)
  elseif ky
    ampo = hubbard_2d_ky(; Nx=Nx, Ny=Ny, t=t, U=U)
  else
    ampo = hubbard_2d(; Nx=Nx, Ny=Ny, yperiodic=yperiodic, t=t, U=U)
  end
  return ampo
end


text:
#2020 Klaus Stricker - Tests for ReadWriteDlm2
#License is MIT: http://julialang.org/license

# rwd2tests_5.jl

# Tests for Tables interface
# ==========================

@testset "5_Stand1 Tables" begin
    mat = [1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"]
    # first, create a MatrixTable from our matrix input
    mattbl = Tables.table(mat)
    # test that the MatrixTable `istable`
    @test Tables.istable(typeof(mattbl))
    # test that it defines row access
    @test Tables.rowaccess(typeof(mattbl))
    @test Tables.rows(mattbl) === mattbl
    # test that it defines column access
    @test Tables.columnaccess(typeof(mattbl))
    @test Tables.columns(mattbl) === mattbl
    # test that we can access the first "column" of our matrix table by column name
    @test mattbl.Column1 == [1,2,3]
    # test our `Tables.AbstractColumns` interface methods
    @test Tables.getcolumn(mattbl, :Column1) == [1,2,3]
    @test Tables.getcolumn(mattbl, 1) == [1,2,3]
    @test Tables.columnnames(mattbl) == [:Column1, :Column2, :Column3]
    # now let's iterate our MatrixTable to get our first MatrixRow
    matrow = first(mattbl)
    @test eltype(mattbl) == typeof(matrow)
    # now we can test our `Tables.AbstractRow` interface methods on our MatrixRow
    @test matrow.Column1 == 1
    @test Tables.getcolumn(matrow, :Column1) == 1
    @test Tables.getcolumn(matrow, 1) == 1
    @test propertynames(mattbl) == propertynames(matrow) == [:Column1, :Column2, :Column3]
end

@testset "5_Stand2 Tables" begin
    rt = [(a=1, b=4.0, c="7"), (a=2, b=5.0, c="8"), (a=3, b=6.0, c="9")]
    ct = (a=[1,2,3], b=[4.0, 5.0, 6.0])
    # let's turn our row table into a plain Julia Matrix object
    mat = Tables.matrix(rt)
    # test that our matrix came out like we expected
    @test mat[:, 1] == [1, 2, 3]
    @test size(mat) == (3, 3)
    @test eltype(mat) == Any
    # so we successfully consumed a row-oriented table,
    # now let's try with a column-oriented table
    mat2 = Tables.matrix(ct)
    @test eltype(mat2) == Float64
    @test mat2[:, 1] == ct.a

    # now let's take our matrix input, and make a column table out of it
    tbl = Tables.table(mat) |> Tables.columntable
    @test keys(tbl) == (:Column1, :Column2, :Column3)
    @test tbl.Column1 == [1, 2, 3]
    # and same for a row table
    tbl2 = Tables.table(mat2) |> Tables.rowtable
    @test length(tbl2) == 3
    @test map(x->x.Column1, tbl2) == [1.0, 2.0, 3.0]
end

@testset "5_matrix2  Tabl" begin
    mat = Any[1 4.0 true "a"; 2 5.0 false "b"; 3 6.0 true "c"]
    # first, create a MatrixTable from our matrix input
    mattbl = Tables.table(mat)
    # test that the MatrixTable `istable`
    @test Tables.istable(typeof(mattbl))
    # ReadWriteDlm2.matrix2: Matrix{Any}, first row Symbol colnames
    ma2 = ReadWriteDlm2.matrix2(mattbl)
    # test that it defines row access
    @test Tables.rowaccess(typeof(mattbl))
    @test Tables.rows(mattbl) === mattbl
    # test that it defines column access
    @test Tables.columnaccess(typeof(mattbl))
    @test Tables.columns(mattbl) === mattbl
    # test that data fit
    @test mattbl.Column1 == [1, 2, 3]
    @test mattbl.Column2 == [4.0, 5.0, 6.0]
    @test mattbl.Column3 == [true, false, true]
    @test mattbl.Column4 == ["a", "b", "c"]
    @test length(mattbl) === 3
    @test ma2[2:end, 1] == [1, 2, 3]
    @test ma2[2:end, 2] == [4.0, 5.0, 6.0]
    @test ma2[2:end, 3] == [true, false, true]
    @test ma2[2:end, 4] == ["a", "b", "c"]
    @test ma2[1, :] == [:Column1, :Column2, :Column3, :Column4]
    # test our `Tables.AbstractColumns` interface methods
    @test Tables.getcolumn(mattbl, :Column1) == [1, 2, 3]
    @test Tables.getcolumn(mattbl, 1) == [1, 2, 3]
    @test Tables.columnnames(mattbl) == [:Column1, :Column2, :Column3, :Column4]
    # now let's iterate our MatrixTable to get our first MatrixRow
    matrow = first(mattbl)
    @test eltype(mattbl) == typeof(matrow)
    # now we can test our `Tables.AbstractRow` interface methods on our MatrixRow
    @test matrow.Column1 == 1
    @test Tables.getcolumn(matrow, :Column1) == 1
    @test Tables.getcolumn(matrow, 1) == 1
    @test propertynames(mattbl) == propertynames(matrow) == [:Column1, :Column2, :Column3, :Column4]
end

@testset "5_Speci  Tables" begin
    mat = [1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"]
    ct = [Int64, Float64, String]
    vof = ReadWriteDlm2.vecofvec(mat, ct) # take columns from array -> Vector{Vector{ct}}
    # first, create a MatrixTable from our Vector{Vector{T}} input
    mattbl = Tables.table(vof)
    # test that the MatrixTable `istable`
    @test Tables.istable(typeof(mattbl))
    # test that it defines row access
    @test Tables.rowaccess(typeof(mattbl))
    @test Tables.rows(mattbl) === mattbl
    # test that it defines column access
    @test Tables.columnaccess(typeof(mattbl))
    @test Tables.columns(mattbl) === mattbl
    # test that we can access the first "column" of our matrix table by column name
    @test mattbl.Column1 == [1,2,3]
    # test our `Tables.AbstractColumns` interface methods
    @test Tables.getcolumn(mattbl, :Column1) == [1,2,3]
    @test Tables.getcolumn(mattbl, 1) == [1,2,3]
    @test Tables.getcolumn(mattbl, ReadWriteDlm2.MatrixTable, 1, :Column1) == [1,2,3]
    @test Tables.columnnames(mattbl) == [:Column1, :Column2, :Column3]
    @test length(mattbl) === 3
    @test iterate(mattbl)[1] == first(mattbl)
    @test iterate(mattbl)[2] == 2
    @test iterate(mattbl, 4) == nothing
    @test Tables.schema(mattbl) ==
        Tables.Schema([:Column1, :Column2, :Column3], [Int64, Float64, String])
    # now let's iterate our MatrixTable to get our first MatrixRow
    matrow = first(mattbl)
    @test eltype(mattbl) == typeof(matrow)
    # now we can test our `Tables.AbstractRow` interface methods on our MatrixRow
    @test matrow.Column1 == 1
    @test matrow.Column2 == 4.0
    @test Tables.getcolumn(matrow, :Column1) == 1
    @test Tables.getcolumn(matrow, 1) == 1
    @test Tables.getcolumn(matrow, ReadWriteDlm2.MatrixRow, 1, :Column1) == 1
    @test propertynames(mattbl) == propertynames(matrow) == [:Column1, :Column2, :Column3]
end

@testset "5_rwdf   Tables" begin
    ma_in = Any[1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"]
    ct = [Int64, Float64, String]
    # take columns from array -> Vector{Vector{ct}}
    vof = ReadWriteDlm2.vecofvec(ma_in, ct)
    # create a MatrixTable from our Vector{Vector{T}} input
    mattbl = Tables.table(vof)
    @test Tables.istable(typeof(mattbl))
    # take data from Matrixtable mattbl in array
    ma_out = ReadWriteDlm2.mttoarray(mattbl, elt=Any)
    # test input array = output array
    @test ma_out == ma_in

    ma_in = Any[1 4.0 "7"; 2 5.0 "8"; 3 6.0 "9"]
    ct = [Float64, Float64, String]
    # take columns from array -> Vector{Vector{ct}}
    vof = ReadWriteDlm2.vecofvec(ma_in, ct)
    # create a MatrixTable from our Vector{Vector{T}} input
    mattbl = Tables.table(vof)
    @test Tables.istable(typeof(mattbl))
    # take data from Matrixtable mattbl in array
    ma_out = ReadWriteDlm2.mttoarray(mattbl, elt=Union{String, Float64})
    # test input array = output array (different types in column1), but ok
    @test ma_out == ma_in

    cn = [:date, :value_1, :value_2]
    mat = [Date(2017,1,1) 1.4 2;
           Date(2017,1,2) 1.8 3;
           nothing missing 4]
    ct = [Union{Nothing, Date}, Union{Missing, Float64}, Int64]
    cna = reshape(cn, 1, :)
    a = vcat(cna, mat)
    writedlm2("test.csv", a)
    astr = read("test.csv", String)
    cstr = "date;value_1;value_2\n2017-01-01;1,4;2\n2017-01-02;1,8;3\nnothing;na;4\n"
    @test astr == cstr
    rdlm1 = readdlm2("test.csv", tables=true, header=true)
    aro1 = ReadWriteDlm2.mttoarray(rdlm1)
    @test isequal(aro1, mat)
    rdlm2 = readdlm2("test.csv", Union{Missing, Nothing, Date, Float64},
                    tables=true, header=true)
    aro2 = ReadWriteDlm2.mttoarray(rdlm2)
    @test isequal(aro2, mat)
    rm("test.csv")
end

@testset "5_csvdlm Tables" begin
    cn = [:date, :value_1, :value_2]
    mat = [Date(2017,1,1) 1.4 2;
           Date(2017,1,2) 1.8 3]
    ct = [Union{Nothing, Date}, Union{Missing, Float64}, Int64]
    vof = ReadWriteDlm2.vecofvec(mat, ct) # take columns from array -> Vector{Vector{ct}}
    # first, create a MatrixTable from our Vector{Vector{T}} input
    mattdf = Tables.table(vof, header=cn)
    # write CSV
    cna = reshape(Tables.columnnames(mattdf), 1, :)
    amt = ReadWriteDlm2.mttoarray(mattdf)
    a = vcat(cna, amt)
    writedlm2("test1.csv", a)
    writecsv2("test2.csv", a)
    # read CSV / Tables Interface
    @test read("test1.csv", String) ==
    "date;value_1;value_2\n2017-01-01;1,4;2\n2017-01-02;1,8;3\n"
    @test read("test2.csv", String) ==
    "date,value_1,value_2\n2017-01-01,1.4,2\n2017-01-02,1.8,3\n"
    df2input1a = readdlm2("test1.csv", dfheader=true)
    df2input2a = readcsv2("test2.csv", dfheader=true)
    @test string(df2input1a) == string(df2input2a)
    df2input1 = readdlm2("test1.csv", tables=true, header=true)
    df2input2 = readcsv2("test2.csv", tables=true, header=true)
    @test string(df2input1) == string(df2input2)
    @test string(df2input1) == string(df2input1a)
    @test string(df2input2) == string(df2input2a)
end

@testset "5_MisNtg Tables" begin
    cn = [:date, :value_1, :value_2]
    mat = [Date(2017,1,1) 1.4 missing;
           Date(2017,1,2) 1.8 nothing]
    ct = [Date, Float64, Union{Missing, Nothing}]
    vof = ReadWriteDlm2.vecofvec(mat, ct) # take columns from array -> Vector{Vector{ct}}
    # first, create a MatrixTable from our Vector{Vector{T}} input
    mattdf = Tables.table(vof, header=cn)
    # write CSV
    cna = reshape(Tables.columnnames(mattdf), 1, :)
    amt = ReadWriteDlm2.mttoarray(mattdf)
    a = vcat(cna, amt)
    writedlm2("test1.csv", a)
    writecsv2("test2.csv", a)
    # read CSV / Tables Interface
    @test read("test1.csv", String) ==
    "date;value_1;value_2\n2017-01-01;1,4;na\n2017-01-02;1,8;nothing\n"
    @test read("test2.csv", String) ==
    "date,value_1,value_2\n2017-01-01,1.4,na\n2017-01-02,1.8,nothing\n"
    df2input1a = readdlm2("test1.csv", dfheader=true)
    df2input2a = readcsv2("test2.csv", dfheader=true)
    @test string(df2input1a) == string(df2input2a)
    df2input1 = readdlm2("test1.csv", tables=true, header=true)
    df2input2 = readcsv2("test2.csv", tables=true, header=true)
    @test string(df2input1) == string(df2input2)
    @test string(df2input1) == string(df2input1a)
    @test string(df2input2) == string(df2input2a)
end


text:
function POMDPModelTools.render(pomdp::SingleOCPOMDP, step::Union{NamedTuple, Dict})
    scene = state_to_scene(pomdp, step.s)
    AutoViz.render(scene, pomdp.env, [TextOverlay(text=["$(step.a)"], font_size=20)], cam = CarFollowCamera(EGO_ID, 20.0))
end

function animate(t, dt)
    frame_index = Int(floor(t/dt)) + 1
    return render(rec[frame_index-nframes(rec)], env.roadway, [CarFollowingStatsOverlay(1)], cam=cam)
end

function animate_states(pomdp::SingleOCPOMDP, states::Vector{SingleOCState}, overlay=nothing)
    duration = length(states)*pomdp.ΔT
    fps = Int(1/pomdp.ΔT)
    cam = FitToContentCamera(0.)
    function render_states(t, dt)
        stateindex = Int(floor(t/dt)) + 1
        scene = state_to_scene(pomdp, states[stateindex])
        return AutoViz.render(scene, pomdp.env, cam = cam)
    end
    return duration, fps, render_states
end

function animate_hist(pomdp::SingleOCPOMDP, hist::SimHistory, overlays::Vector{SceneOverlay}=SceneOverlay[]) where B
    duration = n_steps(hist)*pomdp.ΔT
    fps = Int(1/pomdp.ΔT)
    cam = FitToContentCamera(0.)
    function render_hist(t, dt)
        stateindex = Int(floor(t/dt)) + 1
        scene = state_to_scene(pomdp, state_hist(hist)[stateindex])
        return AutoViz.render(scene, pomdp.env, overlays, cam = cam)
    end
    return duration, fps, render_hist
end

# function animate_record(rec::SceneRecord, overlay::SceneOverlay)
#     duration =rec.nframes*config.sim_dt
#     fps = Int(1/config.sim_dt)
#     function render_rec(t, dt)
#         frame_index = Int(floor(t/dt)) + 1
#         return render(rec[frame_index-nframes(rec)], env, [overlay], cam=cam)
#     end
#     return duration, fps, render_rec
# end
#
# mutable struct DQNFeatures
#     env::CrosswalkEnv
# end
#
# function AutoViz.render!(rendermodel::RenderModel, env_::DQNFeatures)
#     env = env_.env
#     render!(rendermodel, env.roadway)
#
#     obs = env.obstacles[1]
#     for obs in env.obstacles
#         pts = Array{Float64}(2, obs.npts)
#         for (i, pt) in enumerate(obs.pts)
#             pts[1,i] = pt.x
#             pts[2,i] = pt.y
#         end
#
#         add_instruction!(rendermodel, render_fill_region, (pts, colorant"gray"))
#     end
# end
#
# mutable struct InflateOverlay <: SceneOverlay
# end
#
# function AutoViz.render!(rendermodel::RenderModel, overlay::InflateOverlay, scene::Scene, _env::DQNFeatures)
#     env = _env.env
#     ego = scene[findfirst(EGO_ID, scene)]
#     for veh in scene
#         # increase pedestrian size for rendering
#         if veh.def.class == AgentClass.PEDESTRIAN && is_observable(veh.state, ego.state, env)
#             add_instruction!(rendermodel, render_vehicle, (veh.state.posG.x, veh.state.posG.y, veh.state.posG.θ,
#                                               2,2,COLOR_CAR_OTHER))
#         end
#     end
#
#     front = ego.state.posG + polar(VehicleDef().length/2, ego.state.posG.θ)
#     # check if ego is in front of obstacle
#     if front.x < env.params.roadway_length/2 - env.params.crosswalk_width/2
#         obstacles = ConvexPolygon(4)
#         obstacles.pts = [VecSE2(22, -1.5), VecSE2(28, -1.5), VecSE2(28, -13), VecSE2(22, -13)]
#         obstacles.npts = 4; # a rectangle
#         for obs in obstacles
#             pts = Array{Float64}(2, obstacles.npts)
#             for (i, pt) in enumerate(obstacles.pts)
#                 pts[1,i] = pt.x
#                 pts[2,i] = pt.y
#             end
#
#             add_instruction!(rendermodel, render_fill_region, (pts, colorant"black"))
#         end
#     end
#
#     rendermodel
# end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
@testset "Gamma Correction" begin

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        img = fill(oneunit(T),10,10)
        ret = adjust_histogram(img, GammaCorrection(gamma = 1))
        @test img == ret
        @test eltype(ret) == eltype(img)

        imgp = padarray(img, Fill(zero(T), (2,2)))
        retp = adjust_histogram(imgp, GammaCorrection(gamma = 1))
        @test imgp == retp
        @test eltype(retp) == eltype(imgp)
    end
    # ERROR: MethodError: no method matching ^(::AGray{Normed{UInt8,8}}, ::Float64)
    img = fill(oneunit(AGray{N0f8}),10,10)
    @test_broken adjust_histogram(img, GammaCorrection(gamma = 0.5))

    # ERROR: MethodError: no method matching ^(::AGray{Normed{UInt8,8}}, ::Float64)
    img = fill(oneunit(ARGB{N0f8}),10,10)
    @test_broken adjust_histogram(img, GammaCorrection(gamma = 0.5))

    for T in (RGB{N0f8}, RGB{N0f16}, RGB{Float64})
        img = fill(oneunit(T),10,10)
        ret = adjust_histogram(img, GammaCorrection(gamma = 1))
        @test all(map((i, r) -> isapprox(i, r), img, ret))
        @test eltype(ret) == eltype(img)
        imgp = padarray(img, Fill(zero(T), (2,2)))
        retp = adjust_histogram(imgp, GammaCorrection(gamma = 1))
        @test all(map((i, r) -> isapprox(i, r), imgp, retp))
        @test eltype(retp) == eltype(imgp)
    end

    for T in (Gray{N0f8}, Gray{N0f16}, Gray{Float32}, Gray{Float64})
        img = T.(collect(reshape(1/100:1/100:1, 10, 10)))
        for i = 0.5:0.27:2
            ret = adjust_histogram(img, GammaCorrection(gamma = i))
            @test ret == T.(img .^ i)
        end
    end

    # Since the function returns the same output type as the input type
    # there is an implicit rounding operation when dealing with integer values.
    img = reshape(1:1:100, 10, 10)
    for i = 0.5:0.28:2
        ret = adjust_histogram(img, GammaCorrection(gamma = i))
        @test ret == round.(Int, img .^ i)
    end

end


text:
"""
	writebaytap(data,columnout,fileout,sensor;header,decimal)
Write hourly data in Baytap08 fixed/free file format

**Input**
* data: hourly data (equally sampled) in a dataframe containing columns + timevector in DateTime format.
* columnout: which column of the input data should be exported
* fileout: full name of the output file
* sensor: location of the sensor (longitude <-180,180>, latitude <-90,90>,height, mean gravity in cm/s^2)
* header: string to be written to header file  (by default first data column)

**Example**
```
timeout = collect(DateTime(2010,1,1,0):Dates.Hour(1):DateTime(2010,1,2,3));
gravout = ones(Float64,length(timeout));
gravout[[3,end]] .= NaN;
gravout[[4,7]] .= 10.123;
dataout = DataFrame(datetime=timeout, grav=gravout);
writebaytap(dataout,:grav,
			(14.123,45.888,100.0,982.024), # position+mean gravity
			"test/output/baytap_dataseries.txt",
			header="writebaytap unit test");
```
"""
function writebaytap(data::DataFrame,columnout::Symbol,sensor::Tuple{Float64,Float64,Float64,Float64},
					fileout::String;header="Earth tide station")
	# fixed output precission
	decimal = 2;
	val_per_line = 6;
	# Round data to output precision
	dataout = FileTools.round2output(data,[decimal],[columnout]);
	# get time resolutin (input must be equally sampled)
	time_step = Dates.value(diff(dataout[:datetime][1:2])[1])/3600000;
	open(fileout,"w") do fid
		# Write header
		@printf(fid,"%s\n",header);
		@printf(fid,"%10.5f%10.5f%10.2f%10.4f\n",
				sensor[1],sensor[2],sensor[3],sensor[4]);
		@printf(fid,"%s\n",yyyymmddhh(dataout[:datetime][1]))
		@printf(fid,"%5i%5.4g   9000.D0\n",size(dataout,1),time_step)
		@printf(fid,"(%iF10.%d)\n",val_per_line,decimal)
		# Write data
		c = 1;
		for i = 1:size(dataout,1)
			if isnan(dataout[columnout][i])
				@printf(fid,"   9999.99");
			else
				@printf(fid,"%10.2f",dataout[columnout][i]);
			end
			if c < val_per_line # six values per line
				c += 1;
			else
				@printf(fid,"\n");
				c = 1;
			end
		end
	end
end

"""
auxiliary function to export time string in required format
"""
function yyyymmddhh(timein::DateTime)
	yyyy,mm,dd,hh,mi,ss = datevec(timein)
	return @sprintf("%5i%5i%5i%10.5f",
					yyyy,mm,dd,hh + mi/60 .+ ss/3600.)
end

"""
auxiliary function to convert datetime to date vector
"""
function datevec(timein::DateTime)
	return Dates.value(Dates.Year(timein)),
		   Dates.value(Dates.Month(timein)),
		   Dates.value(Dates.Day(timein)),
		   Dates.value(Dates.Hour(timein)),
		   Dates.value(Dates.Minute(timein)),
		   Dates.value(Dates.Second(timein))
end

"""
	baytap2tsoft(file_results,file_output;file_groups,sitename)
Convert baytap08 results (Earth tides) to TSoft tidal group info (written
parameters can be used in "LOCAT.TSD")

**Input**
* file_results: file containing baytap08 results of gravity tidal analysis
* file_output: output file with tidal groups in tsoft format. Paste these results to TSoft/Tidal parameter set
* file_groups: file containing tidal waves after Tamura (1987) (see /test/input/ folder of this package)
* site: name of the site (identical to TSoft database)
* name: name of the tidal parameters (unique for given site)
* filemode: write new file ("w", default) or append to existing ("a")

**Example**
```
file_results = "test/input/baytap08.out";
file_output = "test/output/baytap2tsoft.txt"
baytap2tsoft(file_results,file_output,site="Cantlay",name="test");
```
"""
function baytap2tsoft(file_results::String,file_output::String;
					  site::String="Site",name::String="Name",
					  file_groups::String=joinpath(dirname(@__DIR__),"test","test","input","baytap_tamura_waves.txt"),
					  filemode::String="w")
	waves = readdlm(file_groups);
	# open file for writting
	# open file for writting + add header
	fo = FileTools.tide2tsoft(file_output,site,name,filemode);
	# write parameters
	open(file_results,"r") do fid
		row = readline(fid);
		while !eof(fid)
			if occursin(">",row)
				if row[2] != 'R'
					wstart,wstop,wfactor,wphase,wname = get_wave_info(row);
					freqinfo = find_freq_info(wstart,wstop,waves);
					@printf(fo,"COMP: %8.6f  %8.6f  %7.5f  %8.4f %s\n",
							freqinfo[1],freqinfo[2],wfactor,wphase,wname);
				end
			end
			row = readline(fid);
		end
	end
	close(fo);
end

"""
Auxiliary function to extract needed results
"""
function get_wave_info(row::String)
	temp = split(row)
	wstart = temp |> x -> split(x[3],"-") |> x -> Meta.parse.(x[1])
	wstop = temp |> x -> split(x[3],"-") |> x -> Meta.parse.(x[2][1:end-1])
	wfactor = Meta.parse(temp[5])
	wphase = Meta.parse(temp[7])
	wname = temp[4];
	return wstart,wstop,wfactor,wphase,wname
end

"""
Auxiliary function to find frequency giving wave number
"""
function find_freq_info(wstart,wstop,waves)
	r1 = findall(waves[:,1].==wstart)
	r2 = findall(waves[:,1].==wstop)
	return (waves[r1[1],3],waves[r2[1],3])
end


text:
using CoordinateTransformations: Translation, LinearMap
using MeshCat
using ColorTypes: RGBA, RGB
using MeshCatMechanisms
using StaticArrays
using LinearAlgebra: BLAS
using RigidBodyDynamics
using PressureFieldContact
using Rotations: RotZ
# using PressureFieldContact.Geometry: eMesh_half_plane, eMesh_box, as_tri_eMesh, as_tet_eMesh
using Test

set_zero_subnormals(true)
BLAS.set_num_threads(1)  # NOTE: comment out this line if using IntelMKL

### Box properties
box_rad = 0.05
box_density = 400.0
c_prop_compliant = ContactProperties(Ē=1.0e6)
i_prop_compliant = InertiaProperties(box_density)
i_prop_rigid     = InertiaProperties(box_density, d=box_rad)
eM_box_rigid     = as_tri_eMesh(eMesh_box(box_rad))
eM_box_compliant = as_tet_eMesh(eMesh_box(box_rad))

mech_scen = MechanismScenario()  # n_quad_rule=2)

### Add planes and boxes
nt_plane  = add_contact!(     mech_scen, "plane", as_tet_eMesh(eMesh_half_plane()),   c_prop=c_prop_compliant)
nt_body_1 = add_body_contact!(mech_scen, "box_1", eM_box_rigid,     i_prop=i_prop_rigid)
nt_body_2 = add_body_contact!(mech_scen, "box_2", eM_box_compliant, i_prop=i_prop_compliant, c_prop=c_prop_compliant)
nt_body_3 = add_body_contact!(mech_scen, "box_3", eM_box_rigid,     i_prop=i_prop_rigid)
nt_body_4 = add_body_contact!(mech_scen, "box_4", eM_box_compliant, i_prop=i_prop_compliant, c_prop=c_prop_compliant)

### Friction
add_friction_regularize!(mech_scen, nt_plane.id,  nt_body_1.id, μd=0.0, χ=2.2, n_quad_rule=2)
add_friction_regularize!(mech_scen, nt_body_1.id, nt_body_2.id, μd=0.2, χ=0.2, n_quad_rule=2)
add_friction_regularize!(mech_scen, nt_body_2.id, nt_body_3.id, μd=0.2, χ=0.2, n_quad_rule=2)
add_friction_regularize!(mech_scen, nt_body_3.id, nt_body_4.id, μd=0.2, χ=0.2, n_quad_rule=2)

finalize!(mech_scen)
set_state_spq!(mech_scen, nt_body_1.joint, trans=SVector(0.0, 0.0,  2*box_rad), w=SVector(0.0, 0.0, 1.0))
set_state_spq!(mech_scen, nt_body_2.joint, trans=SVector(0.0, 0.0,  5*box_rad), w=SVector(0.0, 0.0, 2.0))
set_state_spq!(mech_scen, nt_body_3.joint, trans=SVector(0.0, 0.0,  8*box_rad), w=SVector(0.0, 0.0, 3.0))
set_state_spq!(mech_scen, nt_body_4.joint, trans=SVector(0.0, 0.0, 11*box_rad), w=SVector(0.0, 0.0, 4.0))

if (!@isdefined vis)
    vis = Visualizer();
    (!haskey(ENV, "CI")) && open(vis)
end

### Add meshes to visualizer
color_gray  = RGBA{Float32}(0.5, 0.5, 0.5, 1.0)
color_red   = RGBA{Float32}(1.0, 0.0, 0.0, 1.0)
color_green = RGBA{Float32}(0.0, 1.0, 0.0, 1.0)
color_blue  = RGBA{Float32}(0.0, 0.0, 1.0, 1.0)

mvis = MechanismVisualizer(mech_scen, vis)
set_mesh_visual!(     mvis, mech_scen, nt_plane.id,  color_gray)
set_body_mesh_visual!(mvis, mech_scen, nt_body_1.id, color_red)
set_body_mesh_visual!(mvis, mech_scen, nt_body_2.id, color_blue)
set_body_mesh_visual!(mvis, mech_scen, nt_body_3.id, color_green)
set_body_mesh_visual!(mvis, mech_scen, nt_body_4.id, color_red)
set_configuration!(mech_scen, mvis)

### Run forward dynamics
t_final = 5.0e-0
rr = Radau_for_MechanismScenario(mech_scen)
rr.step.h_max = 0.05
@time data_time, data_state = integrate_scenario_radau(rr, t_final=t_final)
println("Finished compiling and running simulation beginning visualization")

### Move camera
setprop!(vis["/Cameras/default/rotated/<object>"], "zoom", 7)
settransform!(vis["/Cameras/default"], Translation(0.0, 0.0, 0.30) ∘ LinearMap(RotZ(-π * 0.35)))

### Visualization
if !haskey(ENV, "CI")
    # If experimenting with this script, you can delete all lines below the word "Visualization" except for the one below.
    play_recorded_data(mvis, mech_scen, data_time, data_state, slowdown=1.0)
else
    function run_visualizer()
        play_recorded_data(mvis, mech_scen, data_time, data_state, slowdown=1.0)
        return true
    end

    @testset "visualizer" begin
        @test run_visualizer()
    end
end


text:
using FlexPlot
using Test

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


text:
module MakieRecipes

using RecipesBase, AbstractPlotting, MakieLayout
using RecipesBase: @recipe
using AbstractPlotting: Palette, to_color

using RecipesPipeline
using Colors

# ## Palette
# The default palette is defined here.

expand_palette(palette, n = 20; kwargs...) = RGBA.(distinguishable_colors(n, palette; kwargs...))

const wong = copy(AbstractPlotting.wong_colors)
begin
    global wong
    tmp = wong[1]
    wong[1] = wong[2]
    wong[2] = tmp
end
const rwong = expand_palette(wong, 50)
const default_palette = Palette(rwong)

include("bezier.jl")
include("pipeline_integration.jl")
include("attribute_table.jl")
include("recipeplot.jl")

# TODO FIXME
RecipesBase.is_key_supported(::Symbol) = false
AbstractPlotting.plots(la::MakieLayout.LAxis) = plots(la.scene)
# FIXME TODO

function tomakie!(sc::AbstractScene, args...; attrs...)
    RecipesPipeline.recipe_pipeline!(sc, Dict{Symbol, Any}(attrs), args)
end

tomakie!(args...; attrs...) = tomakie!(AbstractPlotting.current_scene(), args...; attrs...)

tomakie(args...; attrs...) = tomakie!(Scene(), args...; attrs...)

export tomakie, tomakie!, recipeplot, recipeplot!

end


text:
"""
	within(obj::AbstractMatrix, groups::Vector{Vector{Vector{Int64}}})

	This function performs the within transformation given a model matrix and
	fixed effects using the method of alternating projections.
"""
function within(obj::AbstractMatrix, D::Vector{Vector{Vector{Int64}}})
	output = mapslices(col -> within(col, D), obj, 1)
	return output
end
function within(obj::AbstractVector, D::Vector{Vector{Vector{Int64}}})
	μ = mean(obj)
	current = copy(obj)
	output = copy(obj)
	er = Inf
	while er > 1e-8
		for dimension ∈ D
			current = copy(output)
			for group ∈ dimension
				output[group] .-= mean(current[group])
			end
		end
		er = norm(output - current)
	end
	output .+= μ
	return output
end

function partialwithin(obj::AbstractMatrix, D::Vector{Vector{Int64}},
	θ::Vector{T}) where T <: Real
	means = mean(obj, 1)
	output = copy(obj)
	for idx ∈ eachindex(θ)
		output[D[idx],:] .-= θ[idx] * mean(obj[D[idx],:], 1)
	end
	return output
end
function partialwithin(obj::AbstractVector, D::Vector{Vector{Int64}},
	θ::Vector{T}) where T <: Real
	output = copy(obj)
	for idx ∈ eachindex(θ)
		output[D[idx]] .-= θ[idx] * mean(obj[D[idx]], 1)
	end
	return output
end


text:
mutable struct Output
    categories::Vector{Symbol}
    quantitative::Vector{Symbol}
    summ_quan::Vector{Function}
    design::Matrix{Float64}
    f::Float64
    feas::Vector{Float64}
    TIF::Vector{Float64}
    elapsed_time::Float64
    neighbourhoods::Vector{Neighbourhood}
    infos::Vector{Vector{String}}
    Output(
        categories,
        quantitative,
        summ_quan,
        design,
        f,
        feas,
        TIF,
        elapsed_time,
        neighbourhoods,
        infos,
    ) = new(
        categories,
        quantitative,
        summ_quan,
        design,
        f,
        feas,
        TIF,
        elapsed_time,
        neighbourhoods,
        infos,
    )
    Output() = new(
        Vector{Vector{Symbol}}(undef, 0),
        Vector{Vector{Symbol}}(undef, 0),
        Function[],
        zeros(Float64, 0, 0),
        zero(Float64),
        Float64[],
        Float64[],
        zero(Float64),
        Neighbourhood[],
        Vector{Vector{String}}(undef, 0),
    )
end