text:
typealias Point Vector{Float64}
typealias DelaunayTesselation SubDivision{Point}

minX = 1.0
maxX = 2.0 - eps(Float64)
deltaX = maxX - minX

"""
    DelaunayTesselation()

Creates an initial triangulation consisting of a maximal triangle and the
corresponding faces.

# Remarks
* The inner face coordinates are the centroid of the triangle.
* The outer face coordinates are located at ``(∞, ∞)``.
"""
function DelaunayTesselation()
  local v1 = convert(Point, [minX, minX])
  local v2 = convert(Point, [maxX, minX])
  local v3 = convert(Point, [minX + 0.5 * deltaX, maxX])
  local innerFace = convert(Point, (1./3.) * (v1 + v2 + v3))
  local outerFace = convert(Point, [Inf, Inf])
  local sd = create_subdivision(v1, v2, v3, innerFace, outerFace)

  return sd
end

"""
    getalltriangles(sd::SubDivision)

Returns all triangles as edge lists.

# Remarks
* The boundary edges are included.
"""
function getalltriangles(sd::DelaunayTesselation)
  local triangles = Vector{Tuple{Int, Int, Int}}()
  local nedges = length(sd.edges)
  for ei in 1:2:nedges
    local e2 = lnext(sd, ei)
    local e3 = lnext(sd, e2)
    if (e2 > ei && e3 > ei)
      push!(triangles, (ei, e2, e3))
    end
  end

  return triangles
end

function getdelaunaycoordinates(sd::DelaunayTesselation)
  local x = Vector{Float64}()
  local y = Vector{Float64}()
  for i in 1:4:length(sd.edges)
    local vo = org(sd, i)
    local vd = dest(sd, i)
    if (vo > 3 && vd > 3)
      push!(x, sd.vertices[vo][1])
      push!(x, sd.vertices[vd][1])
      push!(x, NaN)
      push!(y, sd.vertices[vo][2])
      push!(y, sd.vertices[vd][2])
      push!(y, NaN)
    end
  end

  return x, y
end

"""
    swap!(sd::DelaunayTesselation, ei::Int)

Swaps edge `ei` such that the new edge connects the apexes of the triangles
adjacent to the old edge.

# Remarks
* Specifially, if `a = oprev(ei)` and `b = oprev(sym(ei))` before the swap, then
  the edge `ei` will connect `dest(a)` with `dest(b)` after the swap.
"""
function swap!(sd::DelaunayTesselation, ei::Int)
  local es = sym(ei)
  local a = oprev(sd, ei)
  local b = oprev(sd, es)
  local leftFaceIndex = lface(sd, ei)
  local rightFaceIndex = rface(sd, ei)

  # update vertex cache
  sd.vertexCache[org(sd, ei)] = a
  sd.vertexCache[org(sd, es)] = b

  splice!(sd, ei, a)
  splice!(sd, es, b)
  splice!(sd, ei, (lnext(sd, a)))
  splice!(sd, es, (lnext(sd, b)))
  endpoints!(sd, ei, dest(sd, a), dest(sd, b), leftFaceIndex, rightFaceIndex)

  #update faces of adjacent triangles
  local ln = lnext(sd, ei)
  while (ln != ei)
    sd.edges[sym(rot(ln))].origin = leftFaceIndex
    ln = lnext(sd, ln)
  end
  local lnr = lnext(sd, es)
  while (lnr != es)
    sd.edges[sym(rot(lnr))].origin = rightFaceIndex
    lnr = lnext(sd, lnr)
  end
end

"""
    rightof(sd::DelaunayTesselation, vi::Int, ei::Int)
    rightof(sd::DelaunayTesselation, x::Point, ei::Int)

Computes the location of vertex `vi` relative to line defined by `ei`.

# Remarks
* Internally, the computation is forwarded to
  `GeometricalPredicates.orientation`.

# Returns
* True, if `vi` is located to the right of the edge
"""
function rightof(sd::DelaunayTesselation, x::Point, ei::Int)
  local op = sd.vertices[org(sd, ei)]
  local dp = sd.vertices[dest(sd, ei)]
  local o = GeometricalPredicates.Point2D(op[1], op[2])
  local d = GeometricalPredicates.Point2D(dp[1], dp[2])
  local edge = GeometricalPredicates.Line(o,d)
  local test = GeometricalPredicates.Point2D(x[1], x[2])

  return GeometricalPredicates.orientation(edge, test) == -1
end

function rightof(sd::DelaunayTesselation, vi::Int, ei::Int)
  local x = sd.vertices[vi]
  return rightof(sd, x, ei)
end

function closeto(sd::DelaunayTesselation, vi::Int, vj::Int)
  local vip = sd.vertices[vi]
  local vjp = sd.vertices[vj]
  return vip ≈ vjp
end

function edgevect(sd::DelaunayTesselation, ei::Int)
  local v1 = org(sd, ei)
  local v2 = dest(sd, ei)
  return (sd.vertices[v2])-(sd.vertices[v1])
end

function onlinesegment(
  xi::Point,
  x1::Point,
  x2::Point)
    if (xi ≈ x1) || (xi ≈ x2)
      return true
    end
    local d1i = normalize(xi - x1)
    local di2 = normalize(x2 - xi)
    return d1i ≈ di2
end

function onlinesegment(sd::DelaunayTesselation, vi::Int, v1::Int, v2::Int)
  local xi = sd.vertices[vi]
  local x1 = sd.vertices[v1]
  local x2 = sd.vertices[v2]
  if (vi == 284)
    local d1i = normalize(xi - x1)
    local di2 = normalize(x2 - xi)
  end
  return onlinesegment(xi, x1, x2)
end

function orientation(v1::Point, v2::Point, v3::Point)
  local xa = v1[1]
  local ya = v1[2]
  local xb = v2[1]
  local yb = v2[2]
  local xc = v3[1]
  local yc = v3[2]

  local det = (xa*yb - xa*yc - xb*ya + xb*yc + xc*ya - xc*yb)
  if abs(det) <= 1e-10
    return 0
  elseif det > 0.
    return 1
  else
    return 2
  end
end

function orientation(sd::DelaunayTesselation, v1::Int, v2::Int, v3::Int)
  return orientation(sd.vertices[v1], sd.vertices[v2], sd.vertices[v3])
end

"""
    pointsintersect(p1::Point, q1::Point, p2::Point, q2::Point)

Checks if line segments given by their end points intersect.

# Remarks
* Taken from geeksforgeeks.org/check-if-two-given-line-segments-intersect
"""
function pointsintersect(p1::Point, q1::Point, p2::Point, q2::Point)
  local o1 = orientation(p1,q1,p2)
  local o2 = orientation(p1,q1,q2)
  local o3 = orientation(p2,q2,p1)
  local o4 = orientation(p2,q2,q1)
  if (o1 != o2 && o3 != o4)
    return true
  elseif (o1 == 0 && onlinesegment(p2, p1, q1))
    return true
  elseif (o2 == 0 && onlinesegment(q2, p1, q1))
    return true
  elseif (o3 == 0 && onlinesegment(p1, p2, q2))
    return true
  elseif (o4 == 0 && onlinesegment(q1, p2, q2))
    return true
  else
    return false
  end
end

function pointsintersect(sd::DelaunayTesselation, p1::Int, q1::Int, p2::Int, q2::Int)
  return pointsintersect(sd.vertices[p1],sd.vertices[q1],sd.vertices[p2],sd.vertices[q2])
end

function edgeintersects(sd::DelaunayTesselation, e0::Int, v1::Int, v2::Int)
  return pointsintersect(sd, org(sd, e0), dest(sd, e0), v1, v2)
end

function edgeintersects(sd::DelaunayTesselation, e0::Int, e1::Int)
  local vo = org(sd, e1)
  local vd = dest(sd, e1)
  return edgeintersects( sd, e0, vo, vd)
end

"""
    locatevertex(sd::DelaunayTesselation, x::Point, ei::Int)

Locates point `x` in existing subdivision by traversing the triangulation
starting at edge `ei`.

# Remarks
* The point coordinates must be exactly inside the enclosing triangle.
* The algorithm used here [Brown]_ is more stable than the original version

.. [Brown] Brown, Faigle: A robust efficient algorithm for point location in triangulations
"""
function locatevertex(sd::DelaunayTesselation, x::Point, ei::Int)
  cei = ei
  if rightof(sd, x, ei)
    cei = sym(ei)
  end

  local count = 0

  while true
    count += 1
    if (count > length(sd.edges))
      throw(DomainError())
    end

    if x ≈ sd.vertices[org(sd, cei)] || x ≈ sd.vertices[dest(sd, cei)]
      return cei
    end

    local oni = onext(sd, cei)
    local odp = dprev(sd, cei)
    local isrightofOnext = rightof(sd, x, oni)
    local isrightofDprev = rightof(sd, x, odp)
    local op = convert(Int, !isrightofOnext) + 2 * convert(Int, !isrightofDprev)

    if op == 0
      return cei
    elseif op == 1
      cei = oni
    elseif op == 2
      cei = odp
    else
      local e0 = edgevect(sd, sym(oni))
      local xe = x - (sd.vertices[dest(sd, oni)])
      if (dot(e0,xe)) < 0.
        cei = odp
      else
        cei = oni
      end
    end
  end
end

function insertnewvertex!(sd::DelaunayTesselation, e0::Int, vi::Int)
  # add new edge
  local bi = makeedge!(sd)
  local se = bi
  local oldFace = lface(sd, e0)
  endpoints!(sd, bi, org(sd, e0), vi, oldFace, oldFace)

  # connect new edge to enclosing edges
  splice!(sd, bi, e0)

  local nei = e0
  local bni = bi

  while true
    local bsi = sym(bni)
    bni = connect!(sd, nei, bsi, [0.,0.])
    nei = oprev(sd, bni)
    if (lnext(sd, nei) == se)
      break
    end
  end

  return se, nei
end

"""
    insertpoint!(sd::DelaunayTesselation, x::Point, ei::Int)

Inserts a given point into an existing DelaunayMeshes.

# Remarks
* The algorithm identifies the triangle in the existing triangulation that contains
  the vertex to be inserted, taking <paramref name ="ei"/> as start point for the search.
  It then connects the newly-inserted point to all vertices of the enclosing triangle.
* If an existing vertex is already located at the point coordinates an error
  is thrown.

#Returns
* new edge, last edge of target triangle (lnext(lastEdge) = newEdge)
"""
function insertpoint!(sd::DelaunayTesselation, x::Point, ei::Int)
  # find triangle
  local e0 = locatevertex(sd, x, ei)

  # check if the vertex is already in the triangulation
  local eorg = org(sd, e0)
  local edest = dest(sd, e0)
  if (x ≈ sd.vertices[eorg] || x ≈ sd.vertices[edest])
    error("Vertex $x is already in DelaunayMeshes.")
  else
    # insert point as new vertex
    local vi = length(sd.vertices) + 1
    push!(sd.vertices, x)
    push!(sd.vertexCache, -1)

    # check if the point is on any edge of the triangle
    local onEdgeIndex = Nullable{EdgeIndex}()
    if (onlinesegment(sd, vi, eorg, edest))
      onEdgeIndex = Nullable{EdgeIndex}(e0)
    end
    if (isnull(onEdgeIndex))
      local e1 = lnext(sd, e0)
      if (onlinesegment(sd, vi, org(sd, e1), dest(sd, e1)))
        onEdgeIndex = Nullable{EdgeIndex}(e1)
      else
        local e2 = lnext(sd, e1)
        if (onlinesegment(sd, vi, org(sd, e2), dest(sd, e2)))
          onEdgeIndex = Nullable{EdgeIndex}(e2)
        end
      end
    end
    if (!isnull(onEdgeIndex))
      # if the edge is a segment edge we need to repair segment markers
      local ct = sd.edges[get(onEdgeIndex)].constraintType
      if (ct == NoConstraint)
        local ep = oprev(sd, get(onEdgeIndex))
        deleteedge!(sd, get(onEdgeIndex))
        newEdge = insertnewvertex!(sd, ep, vi)
        return newEdge
      else
        local lp = lprev(sd, get(onEdgeIndex))
        local ln = lnext(sd, get(onEdgeIndex))
        deleteedge!(sd, get(onEdgeIndex))
        newEdge = insertnewvertex!(sd, lp, vi)
        local aln = lnext(sd, lp)
        local apn = lprev(sd, ln)
        setconstraint!(sd, aln, ct)
        setconstraint!(sd, apn, ct)
        return newEdge
      end
    else
      newEdge = insertnewvertex!(sd, e0, vi)
      return newEdge
    end
  end
end

"""
    incircle(sd::DelaunayTesselation, v1::Int, v2::Int, v3::Int, v4::Int)

Tests if the vertex `v4` is strictly inside the circumcircle of the triangle
`v1`-`v2`-`v3`.
"""
function incircle(sd::DelaunayTesselation, v1::Int, v2::Int, v3::Int, v4::Int)
  local x1 = GeometricalPredicates.Point2D(sd.vertices[v1][1], sd.vertices[v1][2])
  local x2 = GeometricalPredicates.Point2D(sd.vertices[v2][1], sd.vertices[v2][2])
  local x3 = GeometricalPredicates.Point2D(sd.vertices[v3][1], sd.vertices[v3][2])
  local x4 = GeometricalPredicates.Point2D(sd.vertices[v4][1], sd.vertices[v4][2])
  local triangle = GeometricalPredicates.Primitive(x1,x2,x3)

  return GeometricalPredicates.incircle(triangle, x4) > 0
end

"""
    restoredelaunay!(sd::DelaunayTesselation, ei::Int, si::Int, xi::Int)

Restores Delaunay constraint for edge structure by swapping edges of adjacent triangles.

# Remarks
* It examines all triangles to check if the new vertex `xi` is inside its
  circumcircle.
* The algorithm cycles through all adjacent edges starting at `startEdge`
  and stopping at `stopEdge`.
* If an edge is swapped, its neighbours are examined as well.
"""
function restoredelaunay!(sd::DelaunayTesselation, startEdge::Int, stopEdge::Int, xi::Int)
  local curEdge = startEdge

  while true
    local t = oprev(sd, curEdge)
    local eo = org(sd, curEdge)
    local ed = dest(sd, curEdge)
    local td = dest(sd, t)

    if (rightof(sd, td, curEdge)
       && incircle(sd, eo, td, ed, xi)
       && !isconstraint(sd, curEdge))
      swap!(sd, curEdge)
      curEdge = oprev(sd, curEdge)
    else
      local on = onext(sd, curEdge)
      if (on == stopEdge)
        break
      end
      curEdge = lprev(sd, on)
    end
  end
end

function push!(sd::DelaunayTesselation, x::Array{Float64}, ei::Int)
    local newVertexIndex = length(sd.vertices) + 1
    stopEi, startEi = insertpoint!(sd, x, ei)
    restoredelaunay!(sd, startEi, stopEi, newVertexIndex)
    #=
    local faultyTriangles = TestHelpers.TestDelaunayness(sd)
    @assert(length(faultyTriangles) == 0)
    =#
    return startEi
end

function push!(sd::DelaunayTesselation, points::Array{Float64, 2})
  local startEdge = 1
  for i in 1:Int(length(points)/2)
    startEdge = push!(sd, [points[i, 1], points[i, 2]], startEdge)
  end
end

# for a given triangle, find edge that intersects with v1-v2
function findintersectionedgeoftriangle(
  sd::DelaunayTesselation, e0::Int, v1::Int, v2::Int)

  local ecur = e0
  while true
    if edgeintersects(sd, ecur, v1, v2)
      return ecur
    end

    ecur = lnext(sd, ecur)

    if (ecur == e0)
      throw(ErrorException("Cannot find intersecting edge for triangle $e0."))
    end
  end
end

"""
    findintersectingedges(sd::DelaunayTesselation, v1::Int, v2::Int, eg::Int)

Finds all edges of the given tesselation that intersect with a virtual edge
from `v1` to `v2`.
"""
function findintersectingedges(sd::DelaunayTesselation, v1::Int, v2::Int, eg::Int)
  # find starting edge
  local e0 = locatevertex(sd, sd.vertices[v1], eg)

  # org should be either e0 or lnext e0
  local es = org(sd, e0) == v1 ? e0 : lnext(sd, e0)

  # we need to find the existing edges that enclose the virtual edge
  local ef = es
  local rightOfTest = rightof(sd, v2, es)
  while true
    local nextEdge = onext(sd, ef)
    local rightOfNext = rightof(sd, v2, nextEdge)
    if rightOfNext && !rightOfTest
      break
    end
    ef = nextEdge
    rightOfTest = rightOfNext
  end

  # if enclosing edge has same dest, virtual edge is already in tesselation
  local intersecting = Vector{Int}()
  local alreadyIn = dest(sd, ef) == v2

  if (alreadyIn)
    push!(intersecting, ef)
  else
    while true
      local en = lnext(sd, ef)
      if dest(sd, en) == v2
        break
      end
      local eintersect = findintersectionedgeoftriangle(sd, en, v1, v2)
      ef = sym(eintersect)
      push!(intersecting, eintersect)
    end
  end

  return intersecting, alreadyIn
end

"""
    checkconvexityquadriliteral(sd::DelaunayTesselation, vertices::Vector{Int})

Checks if the given quadriliteral is strictly convex.

# Remarks
* The end points of the quadriliteral must be given in CCW order
"""
function checkconvexityquadriliteral(
  sd::DelaunayTesselation, vertices::Vector{Int})

  return pointsintersect(sd, vertices[1], vertices[3], vertices[2], vertices[4])
end

"""
    getadjacentquadriliteral(sd::DelaunayTesselation, ei::Int)

Gets vertices of quadriliteral formed by the two triangles adjacent to given edge in ccw order.
"""
function getadjacentquadriliteral(sd::DelaunayTesselation, ei::Int)
  local e1 = lnext(sd, ei)
  local e2 = lnext(sd, e1)
  local v2 = org(sd, ei)
  local v0 = org(sd, e1)
  local v1 = org(sd, e2)
  local es = sym(ei)
  local es1 = lnext(sd, es)
  local v3 = dest(sd, es1)
  [v0, v1, v2, v3]
end

"""
    removeintersectingedges!(
      sd::DelaunayTesselation, v1::Int, v2::Int, intersecting::Vector{Int})

Removes all intersecting edges by iteratively swapping diagonals of enclosing quadriliterals.

# Remarks
* For details see Sloan: A fast algorithm for generating constrained Delaunay triangulations.
* Returns the index of the constraint edge
"""
function removeintersectingedges!(
  sd::DelaunayTesselation,
  v1::Int, v2::Int,
  intersecting0::Vector{Int})

  local intersecting = copy(intersecting0)
  local lastSwapped = Nullable{Int}()

  while (length(intersecting) > 0)
    local i = shift!(intersecting)
    local adjacQuadriVert = getadjacentquadriliteral(sd, i)
    if !checkconvexityquadriliteral(sd, adjacQuadriVert)
      # put it at the back of list
      push!(intersecting, i)
    else
      # swap diagonals
      swap!(sd, i)

      # check if swapped edge intersects with virtual edge and put it on list
      # don't put it on the list if org/dest is either of v1/v2
      local io = org(sd, i)
      local id = dest(sd, i)

      if (io == v1) || (io == v2) || (id == v1) || (id == v2)
        lastSwapped = Nullable{Int}(i)
      else
        if (edgeintersects(sd, i, v1, v2))
          push!(intersecting, i)
        else
          lastSwapped = Nullable{Int}(i)
        end
      end
    end
  end

  return get(lastSwapped)
end

function restoredelaunayOfIntersecting!(
  sd::DelaunayTesselation, intersecting::Vector{Int}, v1::Int, v2::Int)

  while (length(intersecting) > 0)
    local i = shift!(intersecting)

    if ((org(sd, i) != v1) && (dest(sd, i) != v2)) && ((org(sd, i) != v2) && (dest(sd, i) != v1))

      local adj = getadjacentquadriliteral(sd, i)

      # check if one of each vertex is contained in other
      local t1in2 = incircle(sd, adj[1], adj[3], adj[4], adj[2])
      local t2in1 = incircle(sd, adj[1], adj[2], adj[3], adj[4])

      # if so swap and put new edge back on list of intersecting
      if t1in2 || t2in1
        swap!(sd, i)
        push!(intersecting, i)
      end
    end
  end
end

"""
    insertconstraint!(sd::DelaunayTesselation, v1::Int, v2::Int)

Inserts constraint going from `v1` to `v2`.

# Remarks
* The constrained region is assumed to be on the right of the directed
  edge.
"""
function insertconstraint!(sd::DelaunayTesselation, v1::Int, v2::Int)
    intersecting, alreadyIn = findintersectingedges(sd, v1, v2, 1)

    if alreadyIn
      local ci = intersecting[1]
      local co = org(sd, ci)
      local ct = v1 == co ? OnRightSide : OnLeftSide
      setconstraint!(sd, ci, ct)
      return ct == OnRightSide ? ci : sym(ci)
    else
      # edge is not already in - we need to swap
      ci = removeintersectingedges!(sd, v1, v2, intersecting)
      co = org(sd, ci)
      ct = v1 == co ? OnRightSide : OnLeftSide
      setconstraint!(sd, ci, ct)
      restoredelaunayOfIntersecting!(sd, intersecting, v1, v2)
      return ct == OnRightSide ? ci : sym(ci)
    end
end

"""
Represents a triangle. It can be constructed from three vertices or a starting edge.
It finds the edges enclosing the triangle and its inner face.
"""
type Triangle
  vertices::Tuple{VertexIndex, VertexIndex, VertexIndex}
  edges::Tuple{Int, Int, Int}
  face::VertexIndex

  function Triangle(sd::DelaunayTesselation, v1::VertexIndex, v2::VertexIndex, v3::VertexIndex)
    # get starting edge from vertex cache
    local ei = sd.vertexCache[v1]
    local e1 = onext(sd, ei)
    while (e1 != ei && dest(sd, e1) != v2)
      e1 = onext(sd, e1)
    end
    if (dest(sd, e1) != v2)
      error("Vertices $v1 - $v2 - $v3 do not belong to the same triangle")
    end
    local e2 = lnext(sd, e1)
    if (dest(sd, e2) != v3)
      error("Vertices $v1 - $v2 - $v3 do not belong to the same triangle")
    end
    local e3 = lnext(sd, e2)
    local face = dest(sd, rot(e1))
    new((v1, v2, v3), (e1, e2, e3), face)
  end

  function Triangle(sd::DelaunayTesselation, e1::EdgeIndex)
    local e2 = lnext(sd, e1)
    local e3 = lnext(sd, e2)
    local v1 = org(sd, e1)
    local v2 = org(sd, e2)
    local v3 = org(sd, e3)
    local face = dest(sd, rot(e1))
    new((v1, v2, v3), (e1, e2, e3), face)
  end
end

# Iterates all triangles
type TrianglesIteratorState
  curIndex::VertexIndex
  visited::Array{VertexIndex}
  curElement::Nullable{VertexIndex}
end

function findNextIteratorState!(sd::DelaunayTesselation, state::TrianglesIteratorState)
  # are there more triangles?
  state.curElement = Nullable{VertexIndex}()
  state.curIndex = state.curIndex + 2
  while (state.curIndex < length(sd.edges))
    cf = org(sd, state.curIndex)
    if (cf > 0 && cf ∉ state.visited)
      state.curElement = cf
      push!(state.visited, cf)
      break
    end
    state.curIndex = state.curIndex + 2
  end

  return state
end

function Base.start(sd::DelaunayTesselation)
  startState = TrianglesIteratorState(0, Array{VertexIndex, 1}(), Nullable{VertexIndex}())
  state = findNextIteratorState!(sd, startState)

  return state
end

function Base.done(sd::DelaunayTesselation, state::TrianglesIteratorState)
  return isnull(state.curElement)
end

function Base.next(sd::DelaunayTesselation, state::TrianglesIteratorState)
  ci = rot(state.curIndex)
  tri = Triangle(sd, ci)

  findNextIteratorState!(sd, state)

  return tri, state
end


text:
module IntArrays

using Mmap

export
    IntArray,
    IntVector,
    IntMatrix,
    radixsort,
    radixsort!

include("buffer.jl")
include("array.jl")
include("vector.jl")
include("matrix.jl")

end # module

ext: jl


text:
@testset "003" begin
    @test tdpe.e003(4) == 2
    @test tdpe.e003(9) == 3
    @test tdpe.e003(13195) == 29
end

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct CppShapeInferenceResult_HandleShapeAndType <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CppShapeInferenceResult_HandleShapeAndType(; kwargs...)
        obj = new(meta(CppShapeInferenceResult_HandleShapeAndType), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CppShapeInferenceResult_HandleShapeAndType
const __meta_CppShapeInferenceResult_HandleShapeAndType = Ref{ProtoMeta}()
function meta(::Type{CppShapeInferenceResult_HandleShapeAndType})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CppShapeInferenceResult_HandleShapeAndType)
            __meta_CppShapeInferenceResult_HandleShapeAndType[] = target = ProtoMeta(CppShapeInferenceResult_HandleShapeAndType)
            allflds = Pair{Symbol,Union{Type,String}}[:shape => TensorShapeProto, :dtype => Int32, :specialized_type => Int32]
            meta(target, CppShapeInferenceResult_HandleShapeAndType, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CppShapeInferenceResult_HandleShapeAndType[]
    end
end
function Base.getproperty(obj::CppShapeInferenceResult_HandleShapeAndType, name::Symbol)
    if name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    elseif name === :dtype
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :specialized_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct CppShapeInferenceResult_HandleData <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CppShapeInferenceResult_HandleData(; kwargs...)
        obj = new(meta(CppShapeInferenceResult_HandleData), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CppShapeInferenceResult_HandleData
const __meta_CppShapeInferenceResult_HandleData = Ref{ProtoMeta}()
function meta(::Type{CppShapeInferenceResult_HandleData})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CppShapeInferenceResult_HandleData)
            __meta_CppShapeInferenceResult_HandleData[] = target = ProtoMeta(CppShapeInferenceResult_HandleData)
            allflds = Pair{Symbol,Union{Type,String}}[:is_set => Bool, :shape_and_type => Base.Vector{CppShapeInferenceResult_HandleShapeAndType}]
            meta(target, CppShapeInferenceResult_HandleData, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CppShapeInferenceResult_HandleData[]
    end
end
function Base.getproperty(obj::CppShapeInferenceResult_HandleData, name::Symbol)
    if name === :is_set
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :shape_and_type
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{CppShapeInferenceResult_HandleShapeAndType}
    else
        getfield(obj, name)
    end
end

mutable struct CppShapeInferenceResult <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CppShapeInferenceResult(; kwargs...)
        obj = new(meta(CppShapeInferenceResult), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CppShapeInferenceResult
const __meta_CppShapeInferenceResult = Ref{ProtoMeta}()
function meta(::Type{CppShapeInferenceResult})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CppShapeInferenceResult)
            __meta_CppShapeInferenceResult[] = target = ProtoMeta(CppShapeInferenceResult)
            fnum = Int[1,4]
            allflds = Pair{Symbol,Union{Type,String}}[:shape => TensorShapeProto, :handle_data => CppShapeInferenceResult_HandleData]
            meta(target, CppShapeInferenceResult, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CppShapeInferenceResult[]
    end
end
function Base.getproperty(obj::CppShapeInferenceResult, name::Symbol)
    if name === :shape
        return (obj.__protobuf_jl_internal_values[name])::TensorShapeProto
    elseif name === :handle_data
        return (obj.__protobuf_jl_internal_values[name])::CppShapeInferenceResult_HandleData
    else
        getfield(obj, name)
    end
end

mutable struct CppShapeInferenceInputsNeeded <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function CppShapeInferenceInputsNeeded(; kwargs...)
        obj = new(meta(CppShapeInferenceInputsNeeded), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct CppShapeInferenceInputsNeeded
const __meta_CppShapeInferenceInputsNeeded = Ref{ProtoMeta}()
function meta(::Type{CppShapeInferenceInputsNeeded})
    ProtoBuf.metalock() do
        if !isassigned(__meta_CppShapeInferenceInputsNeeded)
            __meta_CppShapeInferenceInputsNeeded[] = target = ProtoMeta(CppShapeInferenceInputsNeeded)
            pack = Symbol[:input_tensors_needed,:input_tensors_as_shapes_needed]
            allflds = Pair{Symbol,Union{Type,String}}[:input_tensors_needed => Base.Vector{Int32}, :input_tensors_as_shapes_needed => Base.Vector{Int32}]
            meta(target, CppShapeInferenceInputsNeeded, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, pack, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_CppShapeInferenceInputsNeeded[]
    end
end
function Base.getproperty(obj::CppShapeInferenceInputsNeeded, name::Symbol)
    if name === :input_tensors_needed
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int32}
    elseif name === :input_tensors_as_shapes_needed
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Int32}
    else
        getfield(obj, name)
    end
end

export CppShapeInferenceResult_HandleShapeAndType, CppShapeInferenceResult_HandleData, CppShapeInferenceResult, CppShapeInferenceInputsNeeded


text:
"""

# read_stanfit

Rewrite dct to R format in file.

### Method
```julia
par(cmds)
```

### Required arguments
```julia
* `cmds::Array{Base.AbstractCmd,1}`    : Multiple commands to concatenate

or

* `cmd::Base.AbstractCmd`              : Single command to be
* `n::Number`                            inserted n times into cmd

or
* `cmd::Array{String, 1}`              : Array of cmds as Strings
```

"""
function read_stanfit(model::Stanmodel)

  ## Collect the results of a chain in an array ##

  chainarray = Dict[]

  ## Each chain dictionary can contain up to 5 types of results ##

  result_type_files = ["samples", "diagnostics", "optimize", "diagnose", "variational"]
  rtdict = Dict()

  ## tdict contains the arrays of values ##
  tdict = Dict()

  for i in 1:model.nchains
    for res_type in result_type_files
      if isfile("$(model.name)_$(res_type)_$(i).csv")
        ## A result type file for chain i is present ##
        ## Result type diagnose needs special treatment ##
        instream = open("$(model.name)_$(res_type)_$(i).csv")
        if res_type == "diagnose"
          tdict = Dict()
          str = read(instream, String)
          sstr = split(str)
          tdict = merge(tdict, Dict(:stan_major_version => [parse(Int, sstr[4])]))
          tdict = merge(tdict, Dict(:stan_minor_version => [parse(Int, sstr[8])]))
          tdict = merge(tdict, Dict(:stan_patch_version => [parse(Int, sstr[12])]))
          if tdict[:stan_minor_version][1] >= 10
            # Stan minor version >= 10
            sstr_lp = sstr[79]
            sstr_lp = parse(Float64, split(sstr_lp, '=')[2])
            tdict = merge(tdict, Dict(:lp => [sstr_lp]))
            tdict = merge(tdict, Dict(:var_id => [parse(Int, sstr[90])]))
            tdict = merge(tdict, Dict(:value => [parse(Float64, sstr[91])]))
            tdict = merge(tdict, Dict(:model => [parse(Float64, sstr[92])]))
            tdict = merge(tdict, Dict(:finite_dif => [parse(Float64, sstr[93])]))
            tdict = merge(tdict, Dict(:error => [parse(Float64, sstr[94])]))
          else
            instream = open("$(model.name)_$(res_type)_$(i).csv")
            skipchars(isspace, instream, linecomment='#')
            i = 0
            skipchars(isspace, instream, linecomment='#')
            while true
              i += 1
              line = Unicode.normalize(readline(instream), newline2lf=true)
              if i == 1
                tdict = merge(tdict, Dict(:lp => [float(split(line[1:(length(line)-1)], "=")[2])]))
              elseif i == 3
                sa = split(line[1:(length(line)-1)])
                tdict = merge(tdict, Dict(:var_id => [parse(Int, sa[1])], :value => [float(sa[2])]))
                tdict = merge(tdict, Dict(:model => [float(sa[3])], :finite_dif => [float(sa[4])]))
                tdict = merge(tdict, Dict(:error => [float(sa[5])]))
              end
              if eof(instream)
                close(instream)
                break
              end
              skipchars(isspace, instream, linecomment='#')
            end
          end
        else
          #println("Type of result file is $(res_type)")
          tdict = Dict()
          skipchars(isspace, instream, linecomment='#')
          line = Unicode.normalize(readline(instream), newline2lf=true)
          idx = split(strip(line), ",")
          index = [idx[k] for k in 1:length(idx)]
          #res_type == "optimize" && println(index)
          j = 0
          skipchars(isspace, instream, linecomment='#')
          while true
            j += 1
            #skipchars(isspace, instream, linecomment='#')
            line = Unicode.normalize(readline(instream), newline2lf=true)
            flds = Float64[]
            if eof(instream) && length(line) < 2
              #println("EOF detected")
              close(instream)
              #return(tdict)
              break
            else
              #println(split(strip(line)))
              flds = parse.(Float64, (split(strip(line), ",")))
              for k in 1:length(index)
                if j ==1
                  tdict = merge(tdict, Dict(index[k] => [flds[k]]))
                else
                  tdict[index[k]] = push!(tdict[index[k]], flds[k])
                end
              end
            end
          end
        end

      end

      ## End of processing result type file ##
      ## If any keys were found, merge it in the rtdict ##

      if length(keys(tdict)) > 0
        #println("Merging $(res_type) with keys $(keys(tdict))")
        rtdict = merge(rtdict, Dict(res_type => tdict))
        tdict = Dict()
      end
    end

    ## If rtdict has keys, push it to the chain array ##

    if length(keys(rtdict)) > 0
      #println("Pushing the rtdict with keys $(keys(rtdict))")
      push!(chainarray, rtdict)
      rtdict = Dict()
    end
  end
  chainarray
end

#### Create
####   if m.useMamba = true
####     a Mamba::Chains for results or warmup samples
####   else
####     an Array
####   end

function read_stanfit_warmup_samples(m::Stanmodel, diagnostics=false)
  if !m.method.save_warmup
    println("Saving of warmup samples not requested in StanModel. Returning samples.")
  end
  Stan.read_stanfit_samples(m::Stanmodel, diagnostics, true)
end

function read_stanfit_samples(m::Stanmodel, diagnostics=false, warmup_samples=false)

  global a3d, monitors, index, idx, indvec, ftype

  local noofsamples

  ftype = diagnostics ? "diagnostics" : "samples"

  for i in 1:m.nchains
    if isfile("$(m.name)_$(ftype)_$(i).csv")
      noofsamples = 0
      instream = open("$(m.name)_$(ftype)_$(i).csv")
      skipchars(isspace, instream, linecomment='#')
      line = Unicode.normalize(readline(instream), newline2lf=true)
      idx = split(strip(line), ",")
      index = [idx[k] for k in 1:length(idx)]
      if length(m.monitors) == 0
        indvec = 1:length(index)
      else
        indvec = findall(in(m.monitors), index)
      end
      if m.method.save_warmup
        noofsamples = floor(Int, (m.method.num_samples+m.method.num_warmup)/m.method.thin)
      else
        noofsamples = floor(Int, m.method.num_samples/m.method.thin)
      end
      if i == 1
        a3d = fill(0.0, noofsamples, length(indvec), m.nchains)
      end
      #println(size(a3d))
      skipchars(isspace, instream, linecomment='#')
      for j in 1:noofsamples
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = parse.(Float64, (split(strip(line), ",")))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end
    end
  end

  if m.method.save_warmup
    if warmup_samples
      if m.useMamba
        sr = getindex(a3d, [1:1:size(a3d, 1);],
          [1:size(a3d, 2);], [1:size(a3d, 3);])
        Chains(sr, start=1, thin=1, names=idx[indvec],
          chains=[i for i in 1:m.nchains])
      else
        a3d
      end
    else
      skip_warmup_indx = floor(Int, m.method.num_warmup/m.method.thin) + 1
      if m.useMamba
        sr = getindex(a3d, [skip_warmup_indx:1:noofsamples;],
          [1:size(a3d, 2);], [1:size(a3d, 3);])
        Chains(sr, start=skip_warmup_indx, thin=1,
          names=idx[indvec], chains=[i for i in 1:m.nchains])
      else
        a3d
      end
    end
  else
    if m.useMamba
      sr = getindex(a3d, [1:size(a3d, 1);],
        [1:size(a3d, 2);], [1:size(a3d, 3);])
      Chains(sr, start=1, thin=m.mambaThinning, names=idx[indvec],
        chains=[i for i in 1:m.nchains])
    else
      a3d
    end
  end
end

function read_stanfit_variational_samples(m::Stanmodel)

  local a3d, monitors, index, idx, indvec, ftype
  ftype = "variational"

  for i in 1:m.nchains
    if isfile("$(m.name)_$(ftype)_$(i).csv")
      instream = open("$(m.name)_$(ftype)_$(i).csv")
      skipchars(isspace, instream, linecomment='#')
      line = Unicode.normalize(readline(instream), newline2lf=true)
      idx = split(strip(line), ",")
      index = [idx[k] for k in 1:length(idx)]
      if length(m.monitors) == 0
        indvec = 1:length(index)
      else
        indvec = findall(in(m.monitors), index)
      end
      if i == 1
        a3d = fill(0.0, m.method.output_samples, length(indvec), m.nchains)
      end
      skipchars(isspace, instream, linecomment='#')
      for j in 1:m.method.output_samples
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = parse.(Float64, (split(strip(line), ",")))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end
    end
  end

  if m.useMamba
    sr = getindex(a3d, [1:1:size(a3d, 1);], [1:size(a3d, 2);], [1:size(a3d, 3);])
    Chains(sr, start=1, thin=m.mambaThinning, names=idx[indvec],
      chains=[i for i in 1:m.nchains])
  else
    a3d
  end
end


text:
mutable struct TestMutStruct{T}
    x::T
end

const TEST_VERTEX_TYPES_SMALL = [Int8, Int32, Int64, UInt64]
const TEST_VERTEX_TYPES_BIG =
    [ UInt8,
      Int8,
      Int16,
      UInt16,
      Int32,
      UInt32,
      Int64,
      UInt64,
      Int128,
      UInt128,
      BigInt
    ]

const TEST_EDGEVAL_TYPES_SMALL =
    [ Tuple{},
      Tuple{Float32},
      NamedTuple{(:a,), Tuple{Char}},
      NamedTuple{(:fst, :snd), Tuple{Int8, TestMutStruct{Bool}}}
    ]

const TEST_EDGEVAL_TYPES_SINGE_VALUE_SMALL =
    [ Tuple{TestMutStruct{Char}},
      NamedTuple{(:xyz,), Tuple{Tuple{Int8, UInt8}}},
      Tuple{Union{Int8, Float32}}
    ]

const TEST_EDGEVAL_TYPES_NON_SINGE_VALUE_SMALL =
    [ Tuple{},
      NamedTuple{(), Tuple{}},
      NamedTuple{(:one, :two), Tuple{Int, Int}},
      Tuple{Bool, Int16, UInt16}
    ]

#
# TODO needs more different types, also tuples and named tuples, unitful units, Union types, Nothing

const test_zero_edgval_types = [()]

const test_single_edgeval_types =
    [ (Int,),
      (a=Float16, ),
      (Union{Nothing, Char}, ),
      (a=Tuple{Int, Float32}, ),
      (Vector{Int32}, ),
      (a=Vector, ),
      (NamedTuple{(:a, :b), Tuple{Char, UInt16}}, ),
    ]

const test_multi_edgeval_types =
    [ (Int, Float32),
      (a=Float64, b=Float64),
      (Int8, Int16, Int8),
      (a=Vector{Int64}, b=Tuple{Int, Int}, c=Union{Nothing, Int}),
    ]

const test_edgeval_types = union(
        test_zero_edgval_types,
        test_single_edgeval_types,
        test_multi_edgeval_types
    )

function make_testgraphs(G::Type{<:AbstractGraph}; kwargs...)
    return Channel(c -> _make_testgraphs(c, G, kwargs...))
end

# Similar to `rand` but for more types.
# We do not want type piracy, so we need this.
function rand_sample(::Type{<:Rational{T}}) where {T}
    # TODO do we want to avoid division by zero?
    return rand(T) // rand(T)
end

rand_sample(T::Type) = rand(T)

rand_sample(::Type{Nothing}) = nothing
rand_sample(::Type{Missing}) = missing

rand_sample(::Type{TestMutStruct{T}}) where {T} = TestMutStruct(rand_sample(T))

function rand_sample(U::Union)
    i = rand(1:length(propertynames(U)))
    return rand_sample(getproperty(U, propertynames(U)[i]))
end

function rand_sample(T::Type, dims...)
    result = Array{T}(undef, dims)
    for i in eachindex(result)
        result[i] = rand_sample(T)
    end
    return result
end

function rand_sample(T::Type{<:Union{Tuple, NamedTuple}})
    return T(rand_sample(TT) for TT in T.types)
end

function _make_testgraphs(c::Channel, ::Type{SimpleGraph{V}}; kwargs...) where {V}
        g = SimpleGraph{V}(0)
        info = "$(typeof(g)) 0-graph"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(1)
        info = "$(typeof(g)) 1 vert, 0 edges"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(1)
        add_edge!(g, 1, 1)
        info = "$(typeof(g)) 1 vert, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(2)
        info = "$(typeof(g)) 2 verts, 0 edges"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(2)
        add_edge!(g, 1, 2)
        info = "$(typeof(g)) 2 vers, 1 edge"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(2)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 1 edge, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(2)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 1 edge, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(2)
        add_edge!(g, 1, 2)
        add_edge!(g, 1, 1)
        info = "$(typeof(g)) 2 verts, 1 edge, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleGraph{V}(2)
        add_edge!(g, 1, 2)
        add_edge!(g, 1, 1)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 1 edge, 1 loop"
        put!(c, (graph=g, info=info))

        g = complete_graph(V(5))
        info = "$(typeof(g)) complete_graph(5)"
        put!(c, (graph=g, info=info))
end

function _make_testgraphs(c::Channel, ::Type{SimpleDiGraph{V}}; kwargs...) where{V}
        g = SimpleDiGraph{V}(0)
        info = "$(typeof(g)) 0-graph"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(1)
        info = "$(typeof(g)) 1 vert"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(1)
        add_edge!(g, 1, 1)
        info = "$(typeof(g)) 1 vert, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        info = "$(typeof(g)) 2 verts, 0 edge"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 2)
        info = "$(typeof(g)) 2 verts, 1 edge"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 2 edges, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 1)
        add_edge!(g, 1, 2)
        info = "$(typeof(g)) 2 verts, 1 edge, 1 loop"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 1)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 2 edges, 2 loop"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 2)
        add_edge!(g, 2, 1)
        info = "$(typeof(g)) 2 verts, 2 edges"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 2)
        add_edge!(g, 2, 1)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 3 edges, 1 loops"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 2)
        add_edge!(g, 1, 1)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts, 3 edges, 2 loops"
        put!(c, (graph=g, info=info))

        g = SimpleDiGraph{V}(2)
        add_edge!(g, 1, 2)
        add_edge!(g, 2, 1)
        add_edge!(g, 1, 1)
        add_edge!(g, 2, 2)
        info = "$(typeof(g)) 2 verts 4 edges, 1 loop"
        put!(c, (graph=g, info=info))
end

function testgraphA()
    g = EdgeValGraph(CompleteBipartiteGraph(3, 4))
    add_edge!(g, 1, 4, 10.0)
    return g
end

# TODO only works for undirected graphs
# TODO allow for only non-negative or positive weights
function graph_with_randvals(g::AbstractGraph{T}) where {T}
    U = Float64
    resultg = EdgeValGraph{T, U}(nv(g))
    for e in edges(g)
        s, d = src(e), dst(e)
        add_edge!(resultg, s, d, convert(U, randn()))
    end
    return resultg
end

function is_strictly_increasing(list)
    return all(2:length(list)) do i
        @inbounds list[i - 1] < list[i]
    end
end

# TODO these functions are not finished yet
function testset_isvalidgraph(g::EdgeValGraph; undef_edgevalues::Bool=false)
    @testset "Is valid graph" begin
        @testset "g.ne has correct value" begin
            # Each edge, apart from self-loops occurs twice
            should_be_twice = 0
            for (u, list) in enumerate(g.fadjlist)
                for v in list
                    should_be_twice += ifelse(u == v, 2, 1)
                end
            end
            @test 2 * g.ne == should_be_twice
        end

        @testset "Lists have same sizes" begin
            lengths_should_be = length.(g.fadjlist)
            @test all(g.edgevals) do adjlist
                    length.(adjlist) == lengths_should_be
            end
        end

        @testset "Adjacency list are strictly increasing" begin
            @test all(is_strictly_increasing, g.fadjlist)
        end

        @testset "Symmetric Adjanceny list" begin
            issym = true
            for (u, list_u) in enumerate(g.fadjlist)
                for v in list_u
                    @inbounds list_v = g.fadjlist[v]
                    index = searchsortedfirst(list_v, u)
                    issym &= (index <= length(list_v) && list_v[index] == u)
                end
            end
            @test issym
        end

        if !undef_edgevalues && !isempty(g.edgevals)
            @testset "Symmetric values for key $key" for key in keys(g.edgevals)
                value_adjlist = g.edgevals[key]
                issym = true
                for (u, list_u) in enumerate(g.fadjlist)
                    for (index_u, v) in enumerate(list_u)
                        @inbounds list_v = g.fadjlist[v]
                        index_v = searchsortedfirst(list_v, u)
                        @inbounds issym &= value_adjlist[u][index_u] == value_adjlist[v][index_v]
                    end
                end
                @test issym
            end
        end

    end
end

function testset_isvalidgraph(g::EdgeValOutDiGraph; undef_edgevalues::Bool=false)
    @testset "Is valid graph" begin
        @testset "g.ne has correct value" begin
            should_be = mapreduce(length, +, g.fadjlist, init=0)
            @test g.ne == should_be
        end

        @testset "Lists have same sizes" begin
            lengths_should_be = length.(g.fadjlist)
            @test all(g.edgevals) do adjlist
                    length.(adjlist) == lengths_should_be
            end
        end

        @testset "Adjacency list are strictly increasing" begin
            @test all(is_strictly_increasing, g.fadjlist)
        end
    end
end

function testset_isvalidgraph(g::EdgeValDiGraph; undef_edgevalues::Bool=false)
    @testset "Is valid graph" begin
        @testset "g.ne has correct value" begin
            should_be_forwards = mapreduce(length, +, g.fadjlist, init=0)
            should_be_backwards = mapreduce(length, +, g.badjlist, init=0)
            @test g.ne == should_be_forwards == should_be_backwards
        end

        @testset "Lists have same sizes" begin
            lengths_should_be = length.(g.fadjlist)
            @test all(g.edgevals) do adjlist
                    length.(adjlist) == lengths_should_be
            end
            lengths_should_be = length.(g.badjlist)
            @test all(g.redgevals) do adjlist
                    length.(adjlist) == lengths_should_be
            end
        end

        @testset "Adjacency list are strictly increasing" begin
            @test all(is_strictly_increasing, g.fadjlist)
            @test all(is_strictly_increasing, g.badjlist)
        end

        @testset "Symmetric Adjanceny list" begin
            issym = true
            for (u, list_u) in enumerate(g.fadjlist)
                for v in list_u
                    @inbounds list_v = g.badjlist[v]
                    index = searchsortedfirst(list_v, u)
                    issym &= (index <= length(list_v) && list_v[index] == u)
                end
            end
            @test issym
        end

        if !undef_edgevalues && !isempty(g.edgevals)
            @testset "Symmetric values for key $key" for key in keys(g.edgevals)
                value_adjlist = g.edgevals[key]
                value_badjlist = g.redgevals[key]
                issym = true
                for (u, list_u) in enumerate(g.fadjlist)
                    for (index_u, v) in enumerate(list_u)
                        @inbounds list_v = g.badjlist[v]
                        index_v = searchsortedfirst(list_v, u)
                        @inbounds issym &= value_adjlist[u][index_u] == value_badjlist[v][index_v]
                    end
                end
                @test issym
            end
        end
    end
end

function testset_topological_equivalent(g::SimpleGraph, gv::EdgeValGraph)
    @testset "Topological equivalent" begin
        @test g.fadjlist == gv.fadjlist
    end
end

function testset_topological_equivalent(g::SimpleDiGraph, gv::EdgeValOutDiGraph)
    @testset "Topological equivalent" begin
        @test g.fadjlist == gv.fadjlist
    end
end

function testset_topological_equivalent(g::SimpleDiGraph, gv::EdgeValDiGraph)
    @testset "Topological equivalent" begin
        @test g.fadjlist == gv.fadjlist
        @test g.badjlist == gv.badjlist
    end
end

allkeys_for_E_VALS(E_VALS::Type{<:Tuple}) = 1:length(E_VALS.types)
allkeys_for_E_VALS(E_VALS::Type{<:NamedTuple}) = E_VALS.names ∪
                                                 (1:length(E_VALS.names))


text:
using Combinatorics
using MLStyle

function parse_line(l)
    @match match(r"(\w+) would (lose|gain) (\d+) .* (\w+)\.", l).captures begin
        [u, "lose", n, v] => (u, v, -parse(Int, n))
        [u, "gain", n, v] => (u, v,  parse(Int, n))
        _ => nothing
    end
end

const data = parse_line.(readlines())
const g = begin
    gg = Dict{String, Dict{String, Int}}()
    for (u, v, n) in data
        get!(gg, u, Dict())[v] = n
    end
    gg
end

function happiness(names)
    n = length(names)
    sum(1:n) do i
        get(get(g, names[i], Dict()), names[mod1(i+1, n)], 0) + get(get(g, names[i], Dict()), names[mod1(i-1, n)], 0)
    end
end

# Part 1
g |> keys |> collect |> permutations .|> happiness |> maximum |> println

# Part 2
[collect(keys(g)); "Me"] |> permutations .|> happiness |> maximum |> println

max_stars_repo_name: pLOPeGG/JuliaOfCode

max_issues_repo_name: pLOPeGG/JuliaOfCode

max_forks_repo_name: pLOPeGG/JuliaOfCode


text:
# using SymPy

# x, y = @vars x y
# u = x*(1-x)*y*(1-y)
# v = x*(1-x)*y*(1-y)^2

# p = x*(1-x)*y*(1-y)

# ux = diff(u,x)
# uy = diff(u,y)
# vx = diff(v,x)
# vy = diff(v,y)
# px = diff(p,x)
# py = diff(p,y)
# f = u*ux + v*uy - (diff(ux,x)+diff(uy,y)) + px
# g = u*vx + v*vy - (diff(vx,x)+diff(vy,y)) + py
# h = diff(u, x) + diff(v, y)
# println(replace(replace(sympy.julia_code(simplify(f)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(g)), ".^"=>"^"), ".*"=>"*"))
# println(replace(replace(sympy.julia_code(simplify(h)), ".^"=>"^"), ".*"=>"*"))

function u_exact(x,y)
    if y == 0.0
        return 1.0
    else
        return 0.0
    end
end

function v_exact(x,y)
    0.0
end

function p_exact(x,y)
    0.0
end

function ffunc(x, y)
    #-x^2*y*(x - 1)^2*(y - 1)^2*(2*y - 1) + x*y^2*(x - 1)*(2*x - 1)*(y - 1)^2 + x*y*(y - 1) - 2*x*(x - 1) + y*(x - 1)*(y - 1) - 2*y*(y - 1)
    0.0
end
function gfunc(x, y)
    #x^2*y*(x - 1)^2*(y - 1)^3*(3*y - 1) - x*y^2*(x - 1)*(2*x - 1)*(y - 1)^3 + 3*x*y*(x - 1) + 5*x*(x - 1)*(y - 1) + 2*y*(y - 1)^2
    0.0
end

function hfunc(x,y)
    # (y - 1)*(-2*x*y*(x - 1) + x*y - x*(x - 1)*(y - 1) + y*(x - 1))
    0.0
end

using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays
using ADCMEKit

m = 100
n = 100
h = 1/n
nu = placeholder(0.1)

F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc, m, n, h), m, n, h)
F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc, m, n, h), m, n, h)
H = h^2*eval_f_on_fvm_pts(hfunc, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
Laplace = nu * constant(compute_fem_laplace_matrix1(m, n, h))
function compute_residual(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    f5 = -F1
    F = f1 + f2 + f3 + f4 + f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    g5 = -F2
    G = g1 + g2 + g3 + g4 + g5

    H0 = -B * [u;v] + H

    R = [F;G;H0]
    return R
end

function compute_jacobian(S)
    u, v, p = S[1:(m+1)*(n+1)], S[(m+1)*(n+1)+1:2(m+1)*(n+1)], S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    J0 = [Fu Fv
          Gu Gv]
    J = [J0 -B'
        -B spdiag(zeros(size(B,1)))]
end

NT = 5
# S = zeros(m*n+2(m+1)*(n+1), NT+1)
# S[1:m+1, 1] = ones(m+1,1)

bd = bcnode("all", m, n, h)
# bd = [bd; bd .+ (m+1)*(n+1); ((1:m) .+ 2(m+1)*(n+1))]
bd = [bd; bd .+ (m+1)*(n+1); 2*(m+1)*(n+1)+(n-1)*m+1] # only apply Dirichlet to velocity

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    # residual[bd] .= 0.0
    residual = scatter_update(residual, bd, zeros(length(bd)))

    d = J\residual
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

# for i = 1:NT
#     residual = compute_residual(S[:,i])
#     J = compute_jacobian(S[:,i])

#     J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
#     residual[bd] .= 0.0

#     d = J\residual
#     S[:,i+1] = S[:,i] - d
#     @info i, norm(residual)
# end

xy = fem_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
u0 = @. u_exact(x,y)
v0 = @. v_exact(x,y)

xy = fvm_nodes(m, n, h)
x, y = xy[:,1], xy[:,2]
p0 = @. p_exact(x,y)

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, [u0; v0; p0])

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n))

velo = S[end, 1:2*(m+1)*(n+1)]

velo_data = matread("steady_cavity_data.mat")["V"]
sample_size = 20
idx = rand(1:(m+1)*(n+1), sample_size)
idx = [idx; idx .+ (m+1)*(n+1)] # observe velocity in y direction at the same location
observed_data = velo_data[idx]

noise = true
noise_level = 0.01
if noise
    noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
    observed_data = observed_data .* noise_ratio
end

loss = mean((velo[idx] .- observed_data)^2)
loss = loss * 1e6 # use this to improve optimization accuracy

@info run(sess, loss, nu=>0.01)
# lineview()

# max_iter = 100
# sess = Session(); init(sess)
# loss_ = BFGS!(sess, loss, max_iter)
# figure(); semilogy(loss_); title("loss vs. iteration")

# if noise
#     savefig(string("loss_steady_cavityflow",string(n),"_",string(sample_size),"_",string(noise_level),".png"))
# else
#     savefig(string("loss_steady_cavityflow_",string(n),"_",string(sample_size),".png"))
# end


text:
#--- load modules
using .LAP_julia
using Test, Images, Colors
using FileIO
using CSV
using Plots, ImageView: imshow

#--- load random gray ANHIR image
base_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/"
dataset_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/dataset/"
landmark_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/anhir/landmarks/"
loc_table = CSV.read(base_path * "location_table.csv")

train_rows = loc_table[loc_table[:status] .== "training", :]

const TEST_SIZE = 1

random_rows = train_rows[rand(1:size(train_rows, 1), TEST_SIZE), :]

row_used = (random_rows[1, :])[Symbol("Column1")]

source_loc = (random_rows[1, :])[Symbol("Source image")]
target_loc = (random_rows[1, :])[Symbol("Target image")]
source_in = load(dataset_path * source_loc)
target_in = load(dataset_path * target_loc)

gray_source = Gray.(source_in)
gray_target = Gray.(target_in)

#--- load chessboard image1

tile_size = 4
board_size = 8 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]

chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
chessboard = repeat(mini_board, outer=(2, 2))

imshow(chessboard)

# imshow(chessboard);

#---
small_gray_source = imresize(gray_source, ratio = 0.1)

# LAP filters:
filter_size = 21
filter_half_size = (filter_size - 1) / 2
use_one_dim_filters = false
# Calculate separable filters from basis:
sigma = (filter_half_size + 2) / 4
centered_inds = ImageFiltering.centered(-filter_half_size:filter_half_size)
gaus = ImageFiltering.KernelFactors.gaussian((sigma, sigma), (filter_size, filter_size))
gaus = ImageFiltering.KernelFactors.gaussian(sigma, filter_size)
gaus_der = gaus .* centered_inds .* (-1)/sigma^2
gaus_der_flip = reverse(gaus_der, dims=1)

using Plots
# kernel basis:
# 1
kernf1 = kernelfactors((gaus, gaus))
kernp1 = broadcast(*, kernf1...)
surface(-filter_half_size:filter_half_size, -filter_half_size:filter_half_size, kernp1)

# 2
# forward
kernf2 = kernelfactors((gaus, gaus_der_flip))
kernp2 = broadcast(*, kernf2...)
surface(-filter_half_size:filter_half_size, -filter_half_size:filter_half_size, kernp2)
# backward
kernf2 = kernelfactors((gaus, gaus_der))
kernp2 = broadcast(*, kernf2...)
surface(-filter_half_size:filter_half_size, -filter_half_size:filter_half_size, kernp2)

# 3
# forward
kernf3 = kernelfactors((gaus_der_flip, gaus))
kernp3 = broadcast(*, kernf3...)
surface(-filter_half_size:filter_half_size, -filter_half_size:filter_half_size, kernp3)
# backward
kernf3 = kernelfactors((gaus_der, gaus))
kernp3 = broadcast(*, kernf3...)
surface(-filter_half_size:filter_half_size, -filter_half_size:filter_half_size, kernp3)

surface(-filter_half_size:filter_half_size, -filter_half_size:filter_half_size, kernp)

#--- application to images

# kernels:
# 1
@time imgfilt1 = ImageFiltering.imfilter(chessboard, kernf1, "symmetric")
imshow(imgfilt1)
surface(imgfilt1)

# 2
imgfilt2 = ImageFiltering.imfilter(chessboard, kernf2, "symmetric")
imshow(imgfilt2)
surface(imgfilt2)

# 3
imgfilt3 = ImageFiltering.imfilter(chessboard, kernf3, "symmetric")
imshow(imgfilt3)
surface(imgfilt3)

#--- imfilter with prealloc

imgfilt1_prealloc = similar(chessboard)

@time imgfilt1_tmp = ImageFiltering.imfilter!(imgfilt1_prealloc, chessboard, kernf1, "symmetric")
imshow(imgfilt1_prealloc)
surface(imgfilt1_prealloc)

#---

dif = zeros(4, 4, 3)

dif[:, :, 2] = 1:16
dif[:, :, 3] = ones(4, 4).+1

dif

dif = reshape(dif, (:, 3))

dif

#---

A = zeros(2, 2, 3)
b = zeros(2, 3)

A = rand(100, 100, 100)
b = rand(100, 100)

using BenchmarkTools
@benchmark outpin = pinv(A[:, :, 1])*b[:, 1]
@benchmark outqr = A[:, :, 1] \ b[:, 1]
@benchmark outqrpiv = qr(A[:, :, 1], Val(true)) \ b[:, 1]

outpin ≈ outqr

outqrpiv ≈ outpin

A[:, :, 1] = [1 0; 0 1]
A[:, :, 2] = [1 0; 0 1]
A[:, :, 3] = [1 0; 0 1]

b[:, 1] = [1; 1]
b[:, 2] = [0; 2]
b[:, 3] = [0; 4]

A[:, :, 1]\b[:, 1]

save = zeros(100, 100)

@time multi_mat_div!(save, A, b)

save

function multi_mat_div!(save, A, b)
    for k in axes(A)[3]
        @views save[:, k] = A[:, :, k] \ b[:, k]
    end
end

function multi_mat_div(A, b)
    res = zeros(size(A)[1], size(A)[3])
    # Threads.@threads for j in axes(A, 4) check: (https://stackoverflow.com/questions/57678890/batch-matrix-multiplication-in-julia)
    for k in axes(A)[3]
        @views res[:, k] = qr(A[:, :, k], Val(true)) \ b[:, k]
    end
    return transpose(res)
end

@code_warntype multi_mat_div(A, b)

#---

one = ones(12);
two = reshape(one, (3, 4))

one[:] = two

#--- pad with zeros

using ImageFiltering

image_size = size(chessboard)
filter_half_size = 5;

ones(image_size.-(2*filter_half_size))

border_mask = parent(padarray(ones(image_size.-(2*filter_half_size)), Fill(NaN, (filter_half_size, filter_half_size),(filter_half_size, filter_half_size))))

new_board = chessboard .* border_mask

imshow(new_board)

all_coeffs = ones(broadcast((*), image_size...), 3)

new_coeffs = all_coeffs .* reshape(border_mask, (:, 1))

imshow(reshape(new_coeffs[:, 1], image_size))

imshow(reshape(chessboard[border .== 1.0] .= NaN, (170, 170)))

#---

using BenchmarkTools
using LinearAlgebra

A = rand(11, 11);

A[1, :] = ones(11)

k = (-filter_half_size:filter_half_size)

k1 = repeat(k, outer=11)

sum(transpose(A) * k)

dot(A[:], repeat(k, outer=11))

B = A[:]

#---

u_est, coeffs = single_lap(chessboard, chessboard, 3, 15, [15, 15])

#---

ones_arr_1 = ones(window_size[1])
ones_arr_2 = ones(window_size[2])
ones_kernel = kernelfactors((ones_arr_1, ones_arr_2))
# filtering gets a sum of pixels of window size in each coord

pixels = reshape(chessboard, (:, 1))
image_size = size(chessboard)
window_size = [15, 15]

filter_result = similar(chessboard)

ImageFiltering.imfilter!(filter_result, reshape(pixels, image_size), shit, "symmetric")

imshow(filter_result)

#---

include("../src/lap.jl")

lap.window_sum!(filter_result, pixels, image_size, window_size)

imshow(filter_result)

#---

uv = zeros(size(chessboard)..., 2)
uv[:, :, 1] = real(u_est)
uv[:, :, 2] = imag(u_est)

Plots.quiver([1,2,3], [1,2,3], quiver=([NaN,1,1],[NaN,2,3]))

u_est, coeffs = single_lap(chessboard, chessboard, 3, 15, [15, 15])

p1 = figure()
uv_flow = zeros(size(u_est)..., 2)
uv_flow[:, :, 1] = real(u_est)
uv_flow[:, :, 2] = imag(u_est)

skip_count = 20
n = skip_count

trimmed_real = zeros(size(uv_flow[:,:,1]))
trimmed_real[1:n:end, 1:n:end] = uv_flow[1:n:end, 1:n:end, 1]
trimmed_imag = zeros(size(uv_flow[:,:,2]))
trimmed_imag[1:n:end, 1:n:end] = uv_flow[1:n:end, 1:n:end, 2]

PyPlot.quiver(trimmed_real, trimmed_imag)
gcf()

#---
using FileIO: load

board_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/PolyFilter_LAP/test_boards/"

chess_norm = load(board_path * "chessboard.png")
chess_warp = load(board_path * "chess_warp.png")

chess_norm_f = Float32.(chess_norm)
chess_warp_f = Float32.(chess_warp)

u_est, coeffs = single_lap(chess_norm_f, chess_warp_f, 3, 15, [15, 15])

using LinearAlgebra: qr

qr([1 1; 1 1], Val(true)) \ [1; 1]

#--- TEST single_lap
using PyPlot: quiver, gcf, figure

# *************************
# 1) MATLAB RANDOM FLOW
# *************************
board_path = "/Users/MrTrololord/Google_Drive/cvut/bakalarka/PolyFilter_LAP/test_boards/"

# load chessboards from matlab generated flow
chess_norm = load(board_path * "chessboard.png")
chess_warp = load(board_path * "chess_warp.png")

chess_norm_f = Float32.(chess_norm)
chess_warp_f = Float32.(chess_warp)

u_est, coeffs = LAP_julia.registration.single_lap(chess_norm_f, chess_warp_f, 3, 15, [15, 15])

# *************************
# 2) ONE PIXEL IMAGE SHIFT
# *************************

chess_norm = chessboard
chess_norm_f = Float32.(chess_norm)

chess_warp_f = copy(chess_norm_f)
for k in 2:size(chess_warp_f)[1]
    for l in 2:size(chess_warp_f)[2]
        chess_warp_f[k-1, l-1] = chess_warp_f[k, l]
        if l == size(chess_warp_f)[2]
            chess_warp_f[k-1, l] = chess_warp_f[k, l]
        end
        if k == size(chess_warp_f)[1]
            chess_warp_f[k, l-1] = chess_warp_f[k, l]
        end
    end
end

filter_half_size = 2
filter_size = filter_half_size * 2 + 1

u_est, coeffs = LAP_julia.registration.single_lap(chess_norm_f, chess_warp_f, 3, filter_half_size, [filter_size, filter_size])

# using ImageFiltering: Pad
#
# @views u_est = u_est[]
# u_est = parent(padarray(u_est, Pad(:replicate, filter_half_size, filter_half_size)))

p1 = figure()
PyPlot.quiver(real(u_est), imag(u_est))
gfc()

#---
using Statistics: mean

average = mean(filter(!isnan, u_est))
sqrt(real(average)^2 + imag(average)^2)

#---
image_size = size(chessboard)

window_size = ones(2) .* filter_size
window_half_size = Int64.((window_size .- 1) ./ 2)

Δ_u = ones(10,10)

middle_vals = Δ_u[window_half_size[1]+1:end-window_half_size[1],
                 window_half_size[2]+1:end-window_half_size[2]]
Δ_u = parent(padarray(middle_vals, Pad(:replicate, window_half_size...)))

Δ_u[3, 3] = NaN

Δ_u

any(isnan.(Δ_u))

#---

u = ones(10,10) .+ 2im .* ones(10,10)
u[3:5, 3:5] .= NaN + 1im .* NaN

u

out = LAP_julia.interpolate.interpolate_nans(u)

mask = .!isnan.(u)

u0 = u
u0[.!mask] .= 0

div_coef = similar(u, Int16)

#--- TEST INPAINT

u = ones(10,10) .+ 2im .* ones(10,10);
u[3:8, 3:8] .= NaN .+ NaN .* 1im;

#mask = .!isnan.(real(u))

out = LAP_julia.inpaint.inpaint_nans!(u)

#---

u = ones(10,10) .+ 2im .* ones(10,10);
u[3:8, 3:8] .= NaN .+ NaN .* 1im;

mask = .!isnan.(real(u))

div_coef = ones(size(u))
using ImageFiltering: kernelfactors, imfilter!, Pad, centered

using ImageFiltering: kernelfactors, imfilter, Pad, centered
mask = ones(10, 10)
mask[3:8, 3:8] .= 0

ker = centered([0.5, 2, 0.5])
kernf = kernelfactors((ker, ker));

out = imfilter(mask, kernf, Pad(:symmetric))

#--- Test Makie

using Makie
using TestImages

img = chessboard

scene = Scene()

img1 = reverse(img, dims=1)
img1 = transpose(img1)

scene = image(img, show_axis = false, interpolate = false)

scene = heatmap(img, colormap = :grayscale)

display(scene)

#--- test PyPlot

tile_size = 4
board_size = 8 # must be even
mini_board = [zeros(tile_size, tile_size) ones(tile_size, tile_size);
              ones(tile_size, tile_size) zeros(tile_size, tile_size)]

chessboard = repeat(mini_board, outer=(convert(Integer,(board_size/2)), convert(Integer,(board_size/2))))
chessboard = repeat(mini_board, outer=(2, 2))

img  = chessboard

using Plots

plotly()

Plots.heatmap(img; hover=img)
gui()

###
using PlotlyJS

PlotlyJS.plot(heatmap(z = img))

using Plots
plotly()
data = 1:5
plt = plot(data; hover=data )
display(plt)


text:
using StrideArraysCore, ThreadingUtilities, Aqua
using Test

function closeopensum(x)
    s = zero(eltype(x))
    @inbounds @simd for i ∈ StrideArraysCore.CloseOpen(length(x))
        s += x[i+1]
    end
    s
end
function closeopensumfastmath(x)
    s = zero(eltype(x))
    @inbounds @fastmath for i ∈ StrideArraysCore.CloseOpen(length(x))
        s += x[i+1]
    end
    s
end

@testset "StrideArraysCore.jl" begin

    Aqua.test_all(StrideArraysCore)

    @testset "StrideArrays Basic" begin
        A = rand(100, 100);
        B = copy(A);
        C = StrideArraysCore.PtrArray(A);
        GC.@preserve A begin
            @test closeopensum(C) == closeopensum(A)
            @test closeopensumfastmath(C) == closeopensumfastmath(A)
            @test sum(C) == sum(A)
            @test closeopensum(C) ≈ closeopensumfastmath(C) ≈ sum(C)
            @test A == B
            C .*= 3;
            @test A == 3 .* B
            @test C == A
            @test C == 3 .* B

            D = copy(A);
            Cslice = view(C, 23:48, 17:89)
            Cslice .= 2;
            @test D != C
            D[23:48,17:89] .= 2;
            @test D == C
            @test C === view(C, :, :)
            @test @inferred(size(view(C, StaticInt(1):StaticInt(8), :), 1)) === 8
            @test @inferred(size(view(C, StaticInt(1):StaticInt(8), :), StaticInt(1))) === StaticInt(8)
            @test @inferred(StrideArraysCore.size(view(C, StaticInt(1):StaticInt(8), :), 1)) === 8
            @test @inferred(StrideArraysCore.size(view(C, StaticInt(1):StaticInt(8), :), StaticInt(1))) === StaticInt(8)

            @test C  isa PtrArray
            @test C' isa PtrArray
            @test @inferred(permutedims(C, Val((2,1)))) isa PtrArray
            @test @inferred(adjoint(C)) == D'
            @test @inferred(axes(C)) == axes(D)
            @test @inferred(axes(C')) == axes(D')
            @test @inferred(eachindex(view(C, :, 2:6))) == 1:(5*size(C,1))
            @test @inferred(eachindex(view(C', 2:6, :)')) == 1:(5*size(C,1))
            @test @inferred(eachindex(view(C, 2:6, :))) == CartesianIndices((5, size(C,2)))
            @test @inferred(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), StaticInt(1))) === StrideArraysCore.CloseOpen(StaticInt(5))
            @test @inferred(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), StaticInt(1))) === StrideArraysCore.CloseOpen(StaticInt(5))
            @test @inferred(length(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), StaticInt(1)))) === StaticInt(5)
            @test @inferred(length(axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :), 1))) === 5
            ax1, ax2 = axes(StrideArraysCore.zview(C, StaticInt(2):StaticInt(6), :))
            @test StrideArraysCore.ArrayInterface.known_length(ax1) == 5
            @test StrideArraysCore.ArrayInterface.known_length(ax2) === nothing
            @test StrideArraysCore.ArrayInterface.known_first(ax1) == 0
            @test StrideArraysCore.ArrayInterface.known_first(ax2) == 0
            @test StrideArraysCore.ArrayInterface.known_step(ax1) == 1
            @test StrideArraysCore.ArrayInterface.known_step(ax2) == 1
            @test StrideArraysCore.ArrayInterface.known_last(ax1) == 4
            @test StrideArraysCore.ArrayInterface.known_last(ax2) === nothing
        end
        W = rand(2,3,4);
        X = PtrArray(W);
        GC.@preserve W begin
            @test W == X
            @test permutedims(W, (1,2,3)) == @inferred(permutedims(X, Val((1,2,3))))
            @test permutedims(W, (1,3,2)) == @inferred(permutedims(X, Val((1,3,2))))
            @test permutedims(W, (2,1,3)) == @inferred(permutedims(X, Val((2,1,3))))
            @test permutedims(W, (2,3,1)) == @inferred(permutedims(X, Val((2,3,1))))
            @test permutedims(W, (3,1,2)) == @inferred(permutedims(X, Val((3,1,2))))
            @test permutedims(W, (3,2,1)) == @inferred(permutedims(X, Val((3,2,1))))
            @test_throws BoundsError X[length(X) + 1]
            @test_throws BoundsError X[-4]
            @test_throws BoundsError X[2,5,3]
        end
        @test X === PtrArray(pointer(X), size(X))
        y = rand(77);
        GC.@preserve y begin
            z = PtrArray(y);
            @test y == z
            @test y' == z'
            @test PtrArray(y') === z'
            @test PtrArray(transpose(y)) === transpose(PtrArray(y))
            @test pointer(y) === pointer(z)
            @test_throws BoundsError z[-8]
            @test_throws BoundsError z[88]
        end
    end
    @testset "ThreadingUtilities" begin
        xu = zeros(UInt, 100);
        x = rand(100); y = rand(100); z = rand(100);
        t = (x,y,z)
        pt, gt = StrideArraysCore.object_and_preserve(t)
        greet = ["Hello", "world", "hang", "in", "there"]
        pg, gg = StrideArraysCore.object_and_preserve(greet)
        ph, gh = StrideArraysCore.object_and_preserve(greet[1])
        pht, ght = StrideArraysCore.object_and_preserve((3,greet[1]))
        GC.@preserve xu gt gg gh ght begin
            ThreadingUtilities.store!(pointer(xu), pt, 0)
            @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(pt), 0)) === (sizeof(UInt), t)
            offset = sizeof(UInt)
            for a ∈ t
                _p, g = StrideArraysCore.object_and_preserve(a)
                ThreadingUtilities.store!(pointer(xu), _p, offset)
                offset, p = @inferred(ThreadingUtilities.load(pointer(xu), typeof(_p), offset))
                @test g === a
                @test p isa PtrArray
                @test pointer(p) === pointer(a)
                @test p !== a
                @test p == a
            end
            ThreadingUtilities.store!(pointer(xu), pg, offset)
            @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(pg), offset)) === (offset+sizeof(UInt), greet)
            ThreadingUtilities.store!(pointer(xu), ph, offset)
            @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(ph), offset)) === (offset+sizeof(UInt), greet[1])
            ThreadingUtilities.store!(pointer(xu), pht, offset)
            @test @inferred(ThreadingUtilities.load(pointer(xu), typeof(pht), offset)) === (offset+sizeof(UInt), (3,greet[1]))
        end

        greetsa = StrideArray(greet);
        @test greetsa == greet
        greetsa[1] = "howdy"
        greetsa[5] = "there!"
        @test greet[1] == "howdy"
        @test greet[5] == "there!"
        A = rand(ComplexF64, 4, 5);
        Asa = StrideArray(A);
        @test A == Asa
        Asa[3,4] = 1234.5 + 678.9im;
        @test A[3,4] == 1234.5 + 678.9im;
    end
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct LoadZone <: AggregationTopology
        name::String
        maxactivepower::Float64
        maxreactivepower::Float64
        forecasts::InfrastructureSystems.Forecasts
        internal::InfrastructureSystemsInternal
    end

A collection of buses for electricity price analysis.

# Arguments
- `name::String`
- `maxactivepower::Float64`
- `maxreactivepower::Float64`
- `forecasts::InfrastructureSystems.Forecasts`: internal forecast storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct LoadZone <: AggregationTopology
    name::String
    maxactivepower::Float64
    maxreactivepower::Float64
    "internal forecast storage"
    forecasts::InfrastructureSystems.Forecasts
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function LoadZone(name, maxactivepower, maxreactivepower, forecasts=InfrastructureSystems.Forecasts(), )
    LoadZone(name, maxactivepower, maxreactivepower, forecasts, InfrastructureSystemsInternal(), )
end

function LoadZone(; name, maxactivepower, maxreactivepower, forecasts=InfrastructureSystems.Forecasts(), )
    LoadZone(name, maxactivepower, maxreactivepower, forecasts, )
end

# Constructor for demo purposes; non-functional.
function LoadZone(::Nothing)
    LoadZone(;
        name="init",
        maxactivepower=0.0,
        maxreactivepower=0.0,
        forecasts=InfrastructureSystems.Forecasts(),
    )
end

InfrastructureSystems.get_name(value::LoadZone) = value.name
"""Get LoadZone maxactivepower."""
get_maxactivepower(value::LoadZone) = value.maxactivepower
"""Get LoadZone maxreactivepower."""
get_maxreactivepower(value::LoadZone) = value.maxreactivepower

InfrastructureSystems.get_forecasts(value::LoadZone) = value.forecasts
"""Get LoadZone internal."""
get_internal(value::LoadZone) = value.internal

InfrastructureSystems.set_name!(value::LoadZone, val::String) = value.name = val
"""Set LoadZone maxactivepower."""
set_maxactivepower!(value::LoadZone, val::Float64) = value.maxactivepower = val
"""Set LoadZone maxreactivepower."""
set_maxreactivepower!(value::LoadZone, val::Float64) = value.maxreactivepower = val

InfrastructureSystems.set_forecasts!(value::LoadZone, val::InfrastructureSystems.Forecasts) = value.forecasts = val
"""Set LoadZone internal."""
set_internal!(value::LoadZone, val::InfrastructureSystemsInternal) = value.internal = val


text:
#!/usr/bin/env julia
# -*- coding: utf-8 -*-
# File   : sphero_cylinder_approx.jl
# License: MIT
# Author : Andrei Leonard Nicusan <a.l.nicusan@bham.ac.uk>
# Date   : 21.05.2021

using GLMakie
using Polyhedra
using ColorSchemes
using LinearAlgebra

import CDDLib
lib = CDDLib.Library()

# Return points on a sphere, generated with the method of Markus Deserno from
# https://www.cmu.edu/biolphys/deserno/pdf/sphere_equi.pdf
function half_sphere(;centre = [0. 0. 0.], radius = 1., n = 64)
    a = 4π * radius^2 / 2n
    d = √a

    mθ = floor(π / d)

    dθ = π / mθ
    dϕ = a / dθ

    points = zeros(n, 3)

    # Half sphere without middle circumference
    nc = 0

    for m in 0:mθ / 2 - 2
        θ = π * (m + 0.5) / mθ
        mϕ = floor(2π * sin(θ) / dϕ)
        for p in 0:mϕ - 1
            ϕ = 2π * p / mϕ
            nc += 1

            points[nc, 1] = centre[1] + radius * sin(θ) * cos(ϕ)
            points[nc, 2] = centre[2] + radius * sin(θ) * sin(ϕ)
            points[nc, 3] = centre[3] + radius * cos(θ)
        end
    end

    # Add remaining points on the middle circumference
    inc = 2π / (n - nc)
    for i in 0:n - nc - 1
        nc += 1

        θ = π / 2
        ϕ = i * inc

        points[nc, 1] = centre[1] + radius * sin(θ) * cos(ϕ)
        points[nc, 2] = centre[2] + radius * sin(θ) * sin(ϕ)
        points[nc, 3] = centre[3] + 0.
    end

    points
end

# Construct and return the rotation matrix that rotates a unit vector a onto unit vector b, taken
# from https://math.stackexchange.com/questions/180418
function reorient(a, b)
    v = vec(a) × vec(b)
    vx = [
        0 -v[3] v[2]
        v[3] 0 -v[1]
        -v[2] v[1] 0
    ]

    I(3) + vx + vx * vx ./ (1 + a ⋅ b)
end

function sphero_cylinder(;c1 = [0. -1. -1.], c2 = [4. 4. 0.], r1 = 1., r2 = 1., n = 128)
    hs1 = half_sphere(n = n ÷ 2)
    hs2 = half_sphere(n = n ÷ 2)

    # Re-orient the half-spheres
    v1 = (c1 - c2)
    v1 ./= norm(v1)
    rot1 = reorient([0 0 1], v1)

    for i in 1:size(hs1)[1]
        hs1[i, :] .= rot1 * hs1[i, :]
    end

    v2 = -v1
    rot2 = reorient([0 0 1], v2)

    for i in 1:size(hs2)[1]
        hs2[i, :] .= rot2 * hs2[i, :]
    end

    hs1 = (hs1 .+ c1) .* r1
    hs2 = (hs2 .+ c2) .* r2

    [hs1 ; hs2]
end

# Points on a half-sphere
sphere_points = half_sphere()
sphere_poly = polyhedron(vrep(sphere_points), lib)

# Write the half-sphere vertices to a file
open("half_sphere.csv", "w") do io
    write(io, "Points on a half-sphere polyhedral approximation with 64 vertices\n")
    for i in 1:size(sphere_points)[1]
        for j in 1:size(sphere_points)[2]
            write(io, "$(sphere_points[i, j])\n")
        end
    end
end

# Points on a sphero-cylinder
sc_points = sphero_cylinder()
sc_poly = polyhedron(vrep(sc_points), lib)

# Plot the half-sphere and a sphero-cylinder
cm = ColorSchemes.magma
GLMakie.mesh(Polyhedra.Mesh{3}(sc_poly), color = get(cm, 0.5))
GLMakie.wireframe!(Polyhedra.Mesh{3}(sc_poly))

current_figure()

max_stars_repo_name: anicusan/KonigCell

max_issues_repo_name: anicusan/KonigCell

max_forks_repo_name: anicusan/KonigCell


text:
@testset "StreamKmeans" begin
    seed!(42)
    N = 30
    D = 3
    nInd = 20
    kernel = SqExponentialKernel()
    X = ColVecs(rand(D, N) * 10)
    alg = StreamKmeans(nInd)
    @test repr(alg) == "Streaming Kmeans (m_target=$nInd)"
    test_Zalg(alg; kernel=kernel)
end


text:
## nomodul.jl  Separate Module encapsulation from code
## (c) 2013--2015 David A. van Leeuwen
## Licensed under the MIT software license, see LICENSE.md

## re-load this file during development of the package.

include("modules.jl")

## skip reading types that have been read before
if !isdefined(:TNT)
    include("roctypes.jl")
end

include("includes.jl")


text:
function greet()
    print("Hello World!")
    "Hello World!"
end


text:
module LowLevelJPEG

using AbstractTrees
using ColorTypes
using FixedPointNumbers
using StaticArrays

import Base: ==
export JPEGGray, JPEGGray8, JPEGGray16, JPEGYCbCr, JPEGYCbCr8, JPEGYCbCr16,
       Marker, Markers,
       CodeBits,
       DCSymbol, ACSymbol, DCTCoeff, DCCoeff, ACCoeff,
       DCTBlock,
       QuantizationTable, QuantizationTable8, QuantizationTable16,
       quantize, dequantize,
       HuffmanTreeNode, HuffmanBranchNode,
       HuffmanSymbolNode, DCSymbolNode, ACSymbolNode,
       HuffmanTree, DCHuffmanTree, ACHuffmanTree,
       HuffmanTable, DCHuffmanTable, ACHuffmanTable

u8(io::IO) = read(io, UInt8)
u16(io::IO) = ntoh(read(io, UInt16))
u4u4(io::IO) = (b = u8(io); (b >> 0x4, b & 0xF))

hex(x::Unsigned, pad=2length(x)) = string(x, base=16, pad=pad)

include("colors.jl")
include("markers.jl")
include("codebits.jl")
include("coeffs.jl")
include("blocks.jl")
include("qt.jl")
include("huffman.jl")

end # module


text:
"""
    pop(x::Tuple)

For an `n`-element tuple, returns a tuple `(all_but_last, last)` where `last` is the
final element of the collection and `all_but_last` is a collection containing the first
`n-1` elements.
"""
@inline pop(x::Tuple) = _pop((), x...)
@inline _pop(out::Tuple{}) = error("Can't `pop` on an empty tuple")
@inline _pop(out::Tuple, x) = (out, x)
@inline _pop(out::Tuple, x, y...) = _pop((out..., x), y...)

"""
    pop(x::Tuple, Val{m})

For an `n`-element tuple, returns a tuple `(all_but_last_m, last_m)` where `last_m` is
the final `m` elements of the collection and `all_but_last_m` is a collection containing the
first `n-m` elements.
"""
@generated function pop(x::Tuple, ::Type{Val{N}}) where {N}
    M = length(x.parameters)
    @assert N::Int >= 0
    @assert N <= M

    exprs1 = [:(x[$i]) for i = 1:(M-N)]
    exprs2 = [:(x[$i]) for i = (M-N+1):M]
    quote
        @_inline_meta
        (tuple($(exprs1...)), tuple($(exprs2...)))
    end
end

@inline function same_indices(iter)
    #=i1 = indices(first(iter))
    for i ∈ iter
        if i1 != indices(i)
            return false
        end
    end
    return true =#
#=
    # Optimize
    s = start(iter)
    if done(iter, s)
        return true
    end
    (x, s) = next(iter, s)
    i1 = indices(x)
    while !done(iter, s)
        (x, s) = next(iter,s)
        if _same_indices(indices(x), i1)
            continue
        else
            return false
        end
    end =#
    return true
end

@inline _same_indices(inds1::Tuple{Any}, inds2::Tuple{Any}) = _same_index(inds1[1], inds2[1])
@inline function _same_indices(inds1::Tuple{Any, Any}, inds2::Tuple{Any,Any})
    _same_index(inds1[1], inds2[1]) && _same_index(inds1[2], inds2[2])
end

@inline function _same_index(ind1, ind2)
    if ind1 === ind2
        return true
    end

    i1 = collect(ind1)
    i2 = collect(ind2)

    # They are not equal if they have different indices
    if i1 != i2
        if length(i1) != length(i2)
            for i ∈ i1
                if i ∉ i2
                    return false
                end
            end
        end
    end
    return true
end

@generated function _map(f, x::NTuple{N,Any}) where {N}
    exprs = [:(f(x[$i])) for i = 1:N]
    return quote
        @_inline_meta
        tuple($(exprs...))
    end
end
@generated function _map(f, x::NTuple{N,Any}, y::NTuple{N,Any}) where {N}
    exprs = [:(f(x[$i], y[$i])) for i = 1:N]
    return quote
        @_inline_meta
        tuple($(exprs...))
    end
end
@generated function _map(f, x::NTuple{N,Any}, y::NTuple{N,Any}, z::NTuple{N,Any}) where {N}
    exprs = [:(f(x[$i], y[$i], z[$i])) for i = 1:N]
    return quote
        @_propagate_inbounds_meta
        tuple($(exprs...))
    end
end

@pure issingleton(::Type{T}) where {T} = (length(T.types) == 0 && !T.mutable)
@inline issingleton(x) = issingleton(typeof(x))

@inline aresingleton(x::Tuple) = _aresingleton(x...)
@inline _aresingleton(x, y...) = issingleton(x) && _aresingleton(y...)
@inline _aresingleton(x) = issingleton(x)

@inline aretype(::Type{T}, x::Tuple) where {T} = _aretype(T, x...)
@inline _aretype(::Type{T}, x, y...) where {T} = isa(x, T) && _aretype(T, y...)
@inline _aretype(::Type{T}, x) where {T} = isa(x, T)

_firsttype(::Type{Pair{A,B}}) where {A,B} = A
_firsttype(::Pair{A,B}) where {A,B} = A
_secondtype(::Type{Pair{A,B}}) where {A,B} = B
_secondtype(::Pair{A,B}) where {A,B} = B

@pure increment(i::Int) = i + 1
@pure increment(i::Int, j::Int) = i + j
@pure decrement(i::Int) = i - 1
@pure decrement(i::Int, j::Int) = i - j

# Get the tuple element
@inline get_tuple_value(data::Tuple{Vararg{Pair}}, i) = _get_tuple_value(i, data...)
@inline function _get_tuple_value(i, d::Pair)
    if i == d.first
        return d.second
    else
        throw(IndexError("Can't find index $i"))
    end
end
@inline function _get_tuple_value(i, d::Pair, ds::Pair...)
    if i == d.first
        return d.second
    else
        return _get_tuple_value(i, ds...)
    end
end

# Get the tuple pair
@inline get_tuple_pair(data::Tuple{Vararg{Pair}}, i) = _get_tuple_pair(i, data...)
@inline function _get_tuple_pair(i, d::Pair)
    if i == d.first
        return d
    else
        throw(IndexError("Can't find index $i"))
    end
end
@inline function _get_tuple_pair(i, d::Pair, ds::Pair...)
    if i == d.first
        return d
    else
        return _get_tuple_pair(i, ds...)
    end
end


text:
using GeneralizedDynamicsFromData
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra

#### Setup parameters of the ODE, UDE system and training
parameters_one_stable = Float64[1.5, 1.5, 2.0, 2.0]
parameters_two_stable = Float64[3.5, 3.5, 2.0, 2.0]

parameters = parameters_one_stable

u₀ = Float64[2.0, 4.0]
tspan = (0.0f0, 15.0f0)
s = 0.1

loss_weights = Float64[1/2, 1/2]
max_iters = 2000

#### Setup functions for neural network and differential equations
UA = regression_model(
    2, # input_d
    1, # output_d
    16, # hidden_d
    1, # hidden_layers
    tanh, # nonlinearity
    Flux.glorot_normal # initialization
)

real_de!, univ_de = genetic_toggle_switch(parameters, UA)

#### solve the real ODE and create a target time series
real_prob = ODEProblem(real_de!, u₀, tspan)
real_solution = solve(real_prob, Tsit5(), saveat = s)
y = Array(real_solution)

#### set up the UDE problem and training of the neural network
θ₀ = initial_params(UA)
prob_nn = ODEProblem(univ_de, u₀, tspan, θ₀)

first_guess = concrete_solve(prob_nn, Tsit5(), u₀, θ₀, saveat=real_solution.t)
predict(θ) = Array(concrete_solve(prob_nn, Tsit5(), u₀, θ, saveat=real_solution.t))
loss(θ) = polar_loss(θ, y, loss_weights, predict)

#### train the UDE with ADAM
callback = CallbackLog(T=Float64)
opt = Flux.Optimiser(WeightDecay(1e-4), ExpDecay(1e-1,0.4,100,1e-4), ADAM())
res = DiffEqFlux.sciml_train(loss, θ₀, opt, cb=callback, maxiters=max_iters)

trained_parameters = res.minimizer
ŷ = predict(trained_parameters)
ŷ_derivs = UA(ŷ, trained_parameters)

@variables u[1:2]
b = polynomial_basis(u,3)
recp_polys = []
for i = 1:3
    push!(recp_polys, 1/(1+u[1]^i))
    push!(recp_polys, 1/(1+u[2]^i))
end
gb = [b...; recp_polys...]
basis = Basis(gb, u)

opt = SR3(Float64(1e-2), Float64(1e-1))
λ = exp10.(-7:0.1:5)
g(x) = x[1] > 1 ? Inf : norm(x, 2)

Ψ = SINDy(ŷ,
          ŷ_derivs,
          basis,
          λ,
          opt,
          g = g,
          maxiter = 50000,
          normalize = true,
          denoise = true,
          convergence_error = Float64(1e-10))
println(Ψ)
print_equations(Ψ)

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
function FieldHansenThirdBodyLinear(arg0::jint, arg1::jint, arg2::Field)
    return FieldHansenThirdBodyLinear((jint, jint, Field), arg0, arg1, arg2)
end

function compute_init_values(obj::FieldHansenThirdBodyLinear, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return jcall(obj, "computeInitValues", void, (RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function get_derivative(obj::FieldHansenThirdBodyLinear, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "getDerivative", RealFieldElement, (jint, RealFieldElement), arg0, arg1)
end

function get_value(obj::FieldHansenThirdBodyLinear, arg0::jint, arg1::RealFieldElement)
    return jcall(obj, "getValue", RealFieldElement, (jint, RealFieldElement), arg0, arg1)
end


text:
function position_calc(x...)
    _position_calc(Iterators.filter(x->!isa(x, Nothing), x)...)
end

function normal_calc(x::Bool, invert_normals::Bool = false)
    i = invert_normals ? "-" : ""
    if x
        return "$(i)getnormal(position, position_x, position_y, position_z, o_uv);"
    else
        return "vec3(0, 0, $(i)1);"
    end
end

function light_calc(x::Bool)
    if x
        """
        vec3 L      = normalize(o_lightdir);
        vec3 N      = normalize(o_normal);
        vec3 light1 = blinnphong(N, o_camdir, L, color.rgb);
        vec3 light2 = blinnphong(N, o_camdir, -L, color.rgb);
        color       = vec4(ambient * color.rgb + light1 + backlight * light2, color.a);
        """
    else
        ""
    end
end

function _position_calc(
        position_x::MatTypes{T}, position_y::MatTypes{T}, position_z::MatTypes{T}, target::Type{Texture}
    ) where T<:AbstractFloat
    """
    int index1D = index + offseti.x + offseti.y * dims.x + (index/(dims.x-1));
    ivec2 index2D = ind2sub(dims, index1D);
    vec2 index01 = vec2(index2D) / (vec2(dims)-1.0);
    pos = vec3(
        texelFetch(position_x, index2D, 0).x,
        texelFetch(position_y, index2D, 0).x,
        texelFetch(position_z, index2D, 0).x
    );
    """
end

function _position_calc(
        position_x::VectorTypes{T}, position_y::VectorTypes{T}, position_z::MatTypes{T},
        target::Type{Texture}
    ) where T<:AbstractFloat
    """
    int index1D = index + offseti.x + offseti.y * dims.x + (index/(dims.x-1));
    ivec2 index2D = ind2sub(dims, index1D);
    vec2 index01 = vec2(index2D) / (vec2(dims)-1.0);
    pos = vec3(
        texelFetch(position_x, index2D.x, 0).x,
        texelFetch(position_y, index2D.y, 0).x,
        texelFetch(position_z, index2D, 0).x
    );
    """
end

function _position_calc(
        position_xyz::VectorTypes{T}, target::Type{TextureBuffer}
    ) where T <: StaticVector
    "pos = texelFetch(position, index).xyz;"
end

function _position_calc(
        position_xyz::VectorTypes{T}, target::Type{GLBuffer}
    ) where T <: StaticVector
    len = length(T)
    filler = join(ntuple(x->0, 3-len), ", ")
    needs_comma = len != 3 ? ", " : ""
    "pos = vec3(position $needs_comma $filler);"
end

function _position_calc(
        grid::Grid{2}, position_z::MatTypes{T}, target::Type{Texture}
    ) where T<:AbstractFloat
    """
    int index1D = index + offseti.x + offseti.y * dims.x + (index/(dims.x-1));
    ivec2 index2D = ind2sub(dims, index1D);
    vec2 index01 = vec2(index2D) / (vec2(dims)-1.0);
    float height = texture(position_z, index01).x;
    pos = vec3(grid_pos(position, index01), height);
    """
end

@nospecialize
# surface(::Matrix, ::Matrix, ::Matrix)
function draw_surface(main::Tuple{MatTypes{T}, MatTypes{T}, MatTypes{T}}, data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        position_x = main[1] => (Texture, "x position, must be a `Matrix{Float}`")
        position_y = main[2] => (Texture, "y position, must be a `Matrix{Float}`")
        position_z = main[3] => (Texture, "z position, must be a `Matrix{Float}`")
        scale = Vec3f(0) => "scale must be 0, for a surfacemesh"
    end
    return draw_surface(position_z, data)
end

# surface(Vector or Range, Vector or Range, ::Matrix)
function draw_surface(main::Tuple{VectorTypes{T}, VectorTypes{T}, MatTypes{T}}, data::Dict) where T <: AbstractFloat
    @gen_defaults! data begin
        position_x = main[1] => (Texture, "x position, must be a `Vector{Float}`")
        position_y = main[2] => (Texture, "y position, must be a `Vector{Float}`")
        position_z = main[3] => (Texture, "z position, must be a `Matrix{Float}`")
        scale = Vec3f(0) => "scale must be 0, for a surfacemesh"
    end
    return draw_surface(position_z, data)
end

function draw_surface(main, data::Dict)
    primitive = triangle_mesh(Rect2(0f0,0f0,1f0,1f0))
    to_opengl_mesh!(data, primitive)
    @gen_defaults! data begin
        scale = nothing
        position = nothing
        position_x = nothing => Texture
        position_y = nothing => Texture
        position_z = nothing => Texture
        image = nothing => Texture
        shading = true
        normal = shading
        invert_normals = false
        backlight = 0f0
    end
    @gen_defaults! data begin
        color = nothing => Texture
        color_map = nothing => Texture
        color_norm = nothing
        fetch_pixel = false
        matcap = nothing => Texture

        nan_color = RGBAf(1, 0, 0, 1)
        highclip = RGBAf(0, 0, 0, 0)
        lowclip = RGBAf(0, 0, 0, 0)

        uv_scale = Vec2f(1)
        instances = const_lift(x->(size(x,1)-1) * (size(x,2)-1), main) => "number of planes used to render the surface"
        transparency = false
        shader = GLVisualizeShader(
            "fragment_output.frag", "util.vert", "surface.vert",
            "mesh.frag",
            view = Dict(
                "position_calc" => position_calc(position, position_x, position_y, position_z, Texture),
                "normal_calc" => normal_calc(normal, to_value(invert_normals)),
                "light_calc" => light_calc(shading),
                "buffers" => output_buffers(to_value(transparency)),
                "buffer_writes" => output_buffer_writes(to_value(transparency))
            )
        )
    end
    return assemble_shader(data)
end
@specialize


text:
using DiffPrecTest
using Statistics, StatsBase, LinearAlgebra
using SparseArrays
using ProximalBase, CoordinateDescent, CovSel
using Random, Distributions
using JLD

# rep   = parse(Int,ARGS[1])
# ip    = parse(Int,ARGS[2])
# dir = ARGS[3]

pArr = [50, 100, 150]
n = 300
NUM_REP = 1000

for ip=1:4

    p = pArr[ip]
    Random.seed!(54298)

    # generate model
    Ωx = Matrix{Float64}(I, p, p)
    for l=1:p-1
      Ωx[l  , l+1] = 0.6
      Ωx[l+1, l  ] = 0.6
    end
    for l=1:p-2
      Ωx[l  , l+2] = 0.3
      Ωx[l+2, l  ] = 0.3
    end

    # generate Delta
    d = Vector{Float64}(undef, p)
    rand!(Uniform(0.5, 2.5), d)
    d .= sqrt.(d)
    D = Diagonal(d)
    Σ = inv(Symmetric(D * Ωx * D))

    dist_X = MvNormal(convert(Matrix, Σ))

    rp = div(p * (p+1), 2)
    res = Matrix{DiffPrecResultNormal}(undef, rp, NUM_REP)

    @time for rep=1:NUM_REP
        @show p, rep

        # generate data
        Random.seed!(12346 + rep)
        X = rand(dist_X, n)'
        Y = rand(dist_X, n)'

        Sx = Symmetric( X'X / n )
        Sy = Symmetric( Y'Y / n )

        it = 0
        for col=1:p
            for row=col:p
                it = it + 1
                res[it, rep] = DiffPrecTest.estimate(SeparateNormal(), X, Y, row, col; Sx=Sx, Sy=Sy)
            end
        end

    end

    @save "YinXia_res_$(ip).jld" res
end