text:
# integration code for database usage via CloudGraphs.jl

export
  executeQuery,
  listAllVariables,
  getCloudVert,
  getfnctype,
  usecloudgraphsdatalayer!,
  uselocalmemoryonly!,
  standardcloudgraphsetup,
  consoleaskuserfordb,
  registerGeneralVariableTypes!,
  fullLocalGraphCopy!,
  subLocalGraphCopy!,
  removeGenericMarginals!,
  setBackendWorkingSet!,
  setAllDBSolvable!,
  getExVertFromCloud,
  getAllExVertexNeoIDs,
  getPoseExVertexNeoIDs,
  getFirstPose,
  getfirstpose,
  getLastPose,
  copyAllNodes!,
  copyAllEdges!,
  registerCallback!,
  updateFullCloudVertData!,
  #loading frtend generated fg
  getnewvertdict,
  mergeValuesIntoCloudVert!,
  recoverConstraintType,
  populatenewvariablenodes!,
  populatenewfactornodes!,
  updatenewverts!,
  resetentireremotesession,
  appendvertbigdata!,
  # visualization exports
  getPointCloudFromKinect,
  getPointCloudFromBSON,
  deleteServerSession!

"""
    $(SIGNATURES)

Run Neo4j Cypher queries on the cloudGraph database, and return Tuple with the
unparsed (results, loadresponse).
Throws an error if the query fails.
"""
function executeQuery(connection::Neo4j.Connection,
                      query::AS ) where {AS <: AbstractString}
  #
  loadtx = transaction(connection)
  cph = loadtx(query, submit=true)
  if length(cph.errors) > 0 #Uh oh, return legible error
    err = ""
    for er in cph.errors
      err*="-\r\n"
      for (k,v) in er
        err*=" $k: $v\r\n"
      end
    end
    error("Unable to perform Neo4j query:\r\n$err")
  end
  loadresult = commit(loadtx)
  return cph, loadresult
end
executeQuery(cg::CloudGraph, query::AS) where {AS <:AbstractString} = executeQuery(cg.neo4j.connection, query)

"""
    $(SIGNATURES)

Get cloud vertex from Neo4j using CloudGraphs.jl as identified by `session`, `robot`, `user`, and vertex label `vsym::Symbol`.
"""
function getCloudVert(cg::CloudGraph,
                      session::AbstractString,
                      robot::AbstractString,
                      user::AbstractString,
                      vsym::Symbol;
                      bigdata::Bool=false )
  #
  @warn "getCloudVert(cg, sess, sym) will be deprecated, use getCloudVert(cg, sess, sym=sym) instead."
  # query = " and n.solvable=$(ready) and n.label=$(vsym) "
  # query = reqbackendset ? query*" and n.backendset=$(backendset)" : query
  query = "match (n:$(session):$robot:$user) where n.label='$(vsym)' return id(n)"

  cph, = executeQuery(cg, query)
  # loadtx = transaction(cg.neo4j.connection)
  # cph = loadtx(query, submit=true)
  neoid = cph.results[1]["data"][1]["row"][1]
  CloudGraphs.get_vertex(cg, neoid, bigdata)
end

function getCloudVert(cgl::CloudGraph,
        session::AbstractString,
        robot::AbstractString,
        user::AbstractString;
        exvid::NothingUnion{Int}=nothing,
        neoid::NothingUnion{Int}=nothing,
        sym::NothingUnion{Symbol}=nothing,
        bigdata=false  )
  #
  query = "match (n:$(session):$robot:$user) "

  if sym != nothing
    query = query*" where n.label='$(sym)' "
  elseif exvid != nothing
    query = query*" where n.exVertexId=$(exvid) "
  elseif neoid != nothing
    return CloudGraphs.get_vertex(cgl, neoid, bigdata)
  else
    @show sym, neoid,exvid
    error("Cannot list neighbors if no input reference is given")
  end
  query = query*" return id(n)"

  cph, = executeQuery(cgl.neo4j.connection, query)
  neoid = cph.results[1]["data"][1]["row"][1]
  CloudGraphs.get_vertex(cgl, neoid, bigdata)
end

function listAllVariables(cgl::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)
  #
  query = "match (n:$(session):$robot:$user) where not (n:FACTOR) and exists(n.exVertexId) and n.solvable=1 return n.label, n.exVertexId, id(n), labels(n)"
  cph, = executeQuery(cgl.neo4j.connection, query)

  dd = Dict{Symbol, Tuple{Int, Int, Vector{Symbol}}}()
  for metarows in cph.results[1]["data"]
    data = metarows["row"]
    lbls = setdiff(data[4],String[session])
    enl = (data[2], data[3], Symbol.(lbls))
    dd[Symbol(data[1])] = enl
  end
  return dd
end

function getCloudNeighbors(cgl::CloudGraph, cv::CloudVertex; needdata=false, returntype=false)
  neis = get_neighbors(cgl, cv, needdata=false)
end

"""
    $(SIGNATURES)

List neighbors to node in cgl::CloudGraph by returning Dict{Sym}=(exvid, neoid, Symbol[labels]), and can take
any of the three as input node identifier. Not specifying an identifier will result in all Variable nodes
being returned.
"""
function ls(cgl::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString;
      sym::NothingUnion{Symbol}=nothing,
      neoid::NothingUnion{Int64}=nothing,
      exvid::NothingUnion{Int64}=nothing  )
  #

  if sym == nothing && exvid == nothing && neoid == nothing
    # interrupt and just return all variable nodes
    return listAllVariables(cgl, session, robot, user)
  end
  cv = getCloudVert(cgl, session, robot, user, sym=sym, neoid=neoid, exvid=exvid, bigdata=false )

  neis = get_neighbors(cgl, cv, needdata=false)
  dd = Dict{Symbol, Tuple{Int, Int, Vector{Symbol}}}()
  for nei in neis
    lbls = setdiff(nei.labels, String[session])
    dd[Symbol(nei.properties["label"])] = (nei.exVertexId, nei.neo4jNodeId, Symbol.(lbls) )
  end
  return dd
end

# deprecated
function getfnctype(cvl::CloudGraphs.CloudVertex)
  vert = cloudVertex2ExVertex(cvl)
  return getfnctype(vert)
end

function addCloudVert!(fgl::G,
        exvert::Graphs.ExVertex;
        labels::Vector{T}=String[]  ) where {G <: AbstractDFG, T <: AbstractString}
  #
  cv = CloudGraphs.exVertex2CloudVertex(exvert);
  cv.labels = labels
  CloudGraphs.add_vertex!(fgl.cg, cv);
  fgl.cgIDs[exvert.index] = cv.neo4jNodeId
  IncrementalInference.addGraphsVert!(fgl, exvert)
end

# Return Graphs.ExVertex type containing data according to id
function getExVertFromCloud(fgl::G,
                            fgid::Int64;
                            bigdata::Bool=false  ) where G <: AbstractDFG
  #
  neoID = fgl.cgIDs[fgid]
  cvr = CloudGraphs.get_vertex(fgl.cg, neoID, false)
  CloudGraphs.cloudVertex2ExVertex(cvr)
end

function getExVertFromCloud(fgl::G,
                            lbl::Symbol;
                            nt::Symbol=:var,
                            bigdata::Bool=false  ) where G <: AbstractDFG
  # getExVertFromCloud(fgl, fgl.IDs[lbl], bigdata=bigdata)
  getExVertFromCloud(fgl, (nt==:var ? fgl.IDs[lbl] : fgl.fIDs[lbl]), bigdata=bigdata)
end

function updateFullCloudVertData!(fgl::G,
                                  nv::Graphs.ExVertex;
                                  updateMAPest::Bool=false,
                                  bigdata::Bool=false  ) where G <: AbstractDFG
 #
  # TODO -- this get_vertex seems excessive, but we need the CloudVertex

  if(!haskey(fgl.cgIDs, nv.index))
      error("Cannot find $(nv.index) in $(fgl.cgIDs)...")
  end
  neoID = fgl.cgIDs[nv.index]
  # println("updateFullCloudVertData! -- trying to get $(neoID)")
  vert = CloudGraphs.get_vertex(fgl.cg, neoID, false)
  if typeof(getData(nv)) == VariableNodeData && updateMAPest
    mv = getKDEMax(getKDE(nv))
    nv.attributes["MAP_est"] = mv
    # @show nv.attributes["MAP_est"]
  end

  # TODO -- ignoring other properties
  vert.packed = getData(nv) #.attributes["data"]
  for pair in nv.attributes
    if pair[1] != "data"
      vert.properties[pair[1]] = pair[2]
    end
  end

  # also make sure our local copy is updated, need much better refactoring here
  fgl.stateless ? nothing : fgl.g.vertices[nv.index].attributes["data"] = nv.attributes["data"]

  CloudGraphs.update_vertex!(fgl.cg, vert, bigdata)
end

function makeAddCloudEdge!(fgl::G, v1::Graphs.ExVertex, v2::Graphs.ExVertex) where G <: AbstractDFG
  cv1 = CloudGraphs.get_vertex(fgl.cg, fgl.cgIDs[v1.index], false)
  cv2 = CloudGraphs.get_vertex(fgl.cg, fgl.cgIDs[v2.index], false)
  ce = CloudGraphs.CloudEdge(cv1, cv2, "DEPENDENCE");
  retrel = CloudGraphs.add_edge!(fgl.cg, ce);

  # TODO -- keep this edge id in function node data, must refactor
  push!(v2.attributes["data"].edgeIDs, retrel.id) # TODO -- not good way to do this
  updateFullCloudVertData!(fgl, v2)

  IncrementalInference.makeAddEdge!(fgl, v1, v2, saveedgeID=false)
  retrel.id
end

# TODO -- fetching of CloudVertex propably not required, make faster request to @GearsAD
function getCloudOutNeighbors(fgl::G,
                              exVertId::Int64;
                              solvable::Int=1,
                              backendset::Int=1,
                              needdata::Bool=false  ) where G <: AbstractDFG
  #
  # println("Looking for cloud out neighbors")
  cgid = fgl.cgIDs[exVertId]
  cv = CloudGraphs.get_vertex(fgl.cg, cgid, false)
  neighs = CloudGraphs.get_neighbors(fgl.cg, cv)
  neExV = Graphs.ExVertex[]
  for n in neighs
    cgn = CloudGraphs.cloudVertex2ExVertex(n)
    if (
       #cgn.attributes["solvable"] == solvable &&
       #cgn.attributes["backendset"] == backendset &&
       (!needdata || haskey(cgn.attributes, "data") )  )
      push!(neExV, cgn )
    end
  end
  return neExV
end

# return list of neighbors as Graphs.ExVertex type
function getCloudOutNeighbors(fgl::G,
                              vert::Graphs.ExVertex;
                              solvable::Int=1,
                              backendset::Int=1,
                              needdata::Bool=false  ) where G <: AbstractDFG
  # TODO -- test for solvable and backendset here
  getCloudOutNeighbors(fgl, vert.index, solvable=solvable,backendset=backendset, needdata=needdata )
end

function getEdgeFromCloud(fgl::G, id::Int64) where G <: AbstractDFG
  println("getting id=$(id)")
  CloudGraphs.get_edge(fgl.cg, id)
end

function deleteCloudVertex!(fgl::G, vert::Graphs.ExVertex) where G <: AbstractDFG
  neoID = fgl.cgIDs[vert.index]
  cvr = CloudGraphs.get_vertex(fgl.cg, neoID, false)
  CloudGraphs.delete_vertex!(fgl.cg, cvr)
end

function deleteCloudEdge!(fgl::G, edge::CloudEdge) where G <: AbstractDFG
  CloudGraphs.delete_edge!(fgl.cg, edge)
end

function usecloudgraphsdatalayer!()
  IncrementalInference.setdatalayerAPI!(
    addvertex= addCloudVert!,
    getvertex= getExVertFromCloud,
    makeaddedge= makeAddCloudEdge!,
    getedge= getEdgeFromCloud,
    outneighbors= getCloudOutNeighbors,
    updatevertex= updateFullCloudVertData!,
    deletevertex= deleteCloudVertex!,
    deleteedge= deleteCloudEdge!,
    cgEnabled= true )
  nothing
end

function uselocalmemoryonly!()
  IIF.setdatalayerAPI!(
    addvertex= IIF.localapi.addvertex!,
    getvertex= IIF.localapi.getvertex,
    makeaddedge= IIF.localapi.makeaddedge!,
    getedge= IIF.localapi.getedge,
    outneighbors= IIF.localapi.outneighbors,
    updatevertex= IIF.localapi.updatevertex!,
    deletevertex= IIF.localapi.deletevertex!,
    deleteedge= IIF.localapi.deleteedge!,
    cgEnabled= false )
  nothing
end

function getpackedtype(typestring::AS) where {AS <: AbstractString}
  # println("Caesar.getpackedtype($(typestring))")
  eval(Meta.parse(typestring))() # TODO consider caching or better
end

# function should not be necessary, but fixes a minor bug following elimination algorithm
function removeGenericMarginals!(conn)
  loadtx = transaction(conn)
  query = "match (n)-[r]-() where n.packedType = 'IncrementalInference.FunctionNodeData{IncrementalInference.GenericMarginal}' detach delete n,r"
  cph = loadtx(query, submit=true)
  loadresult = commit(loadtx)
  # TODO -- can probably be made better, but should not be necessary in the first place
  loadtx = transaction(conn)
  query = "match (n) where n.packedType = 'IncrementalInference.FunctionNodeData{IncrementalInference.GenericMarginal}' detach delete n"
  cph = loadtx(query, submit=true)
  loadresult = commit(loadtx)
  nothing
end

"""
    $(SIGNATURES)

Get all Neo4j node IDs in current session.
"""
function getAllExVertexNeoIDs(conn::Neo4j.Connection;
        solvable::Int=1,
        backendset::Int=1,

        sessionname::AS="",
        robotname::AS="",
        username::AS="",
        reqbackendset::Bool=true,
        reqSolvable::Bool=true  ) where {AS <: AbstractString}
  #
  sn = length(sessionname) > 0 ? ":"*sessionname : ""
  rn = length(robotname) > 0 ? ":"*robotname : ""
  un = length(username) > 0 ? ":"*username : ""
  query = "match (n$(sn)$(rn)$(un)) where not n:SESSION and not n:MULTISESSION and exists(n.exVertexId)"
  query = reqbackendset || reqSolvable ? query*" and" : query
  query = reqSolvable ? query*" n.solvable=$(solvable)" : query
  query = reqbackendset && reqSolvable ? query*" and" : query
  query = reqbackendset ? query*" n.backendset=$(backendset)" : query
  query = query*" return n.exVertexId, id(n), n.label"

  cph, = executeQuery(conn, query)
  ret = Array{Tuple{Int64,Int64,Symbol},1}()

  # If no results, return empty array...
  if length(cph.results) > 0
      if length(cph.results[1]["data"]) > 0
          @showprogress 1 "Get ExVertex IDs..." for data in cph.results[1]["data"]
            exvid, neoid, sym = data["row"][1], data["row"][2], Symbol(data["row"][3])
            push!(ret, (exvid,neoid,sym)  )
          end
      end
  end

  return ret
end

"""
    $(SIGNATURES)

Build query to fetch sub graph and neighboring nodes.  For example:
FAILS IN SOME CASES
```
match (n0:Hackathon:HexagonalDrive)
where n0.label IN ['x0']
with n0 as m
return m.exVertexId, id(m), m.label
UNION
match (n0:Hackathon:HexagonalDrive)-[:DEPENDENCE]-(n1:HexagonalDrive)
where n0.label IN ['x0']
with n1 as m
return m.exVertexId, id(m), m.label
UNION
match (n0:Hackathon:HexagonalDrive)-[:DEPENDENCE]-(n1:HexagonalDrive)-[:DEPENDENCE]-(n2:HexagonalDrive)
where n0.label IN ['x0']
with n2 as m
return m.exVertexId, id(m), m.label
```
"""
function buildSubGraphIdsQuery(;
            lbls::Vector{AS}=String[""],
            session::AS="",
            robot::AS="",
            user::AS="",
            label::AS="",
            reqSolvable::Bool=true,
            solvable::Int=1,
            reqbackendset::Bool=true,
            backendset::Int=1,
            neighbors::Int=0,
            includeMultisession::Bool=false ) where {AS <: AbstractString}
  #
  sn = length(session) > 0 ? ":"*session : ""
  rn = length(robot) > 0 ? ":"*robot : ""
  un = length(user) > 0 ? ":"*user : ""
  lb = length(label) > 0 ? ":"*label : ""

  query = ""
  for nei in 0:(neighbors)
    query *= "match (n0$(sn)$(rn)$(un)$(lb))"
    outerd = 0
    for d in 1:(nei)
      query *= "-[:DEPENDENCE]-(n$(d)$(sn)$(rn)$(un)$(lb))"
      outerd = d
    end
    query *= " "
    query *= "where n0.label IN ["
    for lbl in lbls
      query *= "'$(lbl)',"
    end
    query = chop(query)*"] "
    if !includeMultisession
      query *= "and not n0:MULTISESSION "
    end
    query *= "with n$(outerd) as m "
    query *= "return m.exVertexId, id(m), m.label"
    query *= nei != (neighbors) ? " UNION " : ""
  end
  return query
end
# qu = buildSubGraphIdsQuery(lbls=["x0";], session="HexagonalDrive", neighbors=3)
# @show qu

"""
    $(SIGNATURES)

Return array of tuples with ExVertex IDs and Neo4j IDs for vertices with label in session.
"""
function getLblExVertexNeoIDs(
        conn::Neo4j.Connection,
        lbls::Vector{AS};
        session::AS="",
        robot::AS="",
        user::AS="",
        label::AS="",
        reqSolvable::Bool=true,
        solvable::Int=1,
        backendset::Int=1,
        reqbackendset::Bool=true,
        neighbors::Int=0,
        includeMultisession::Bool=false) where {AS <: AbstractString}
  #

  query = buildSubGraphIdsQuery(lbls=lbls, session=session, robot=robot, user=user, label=label, neighbors=neighbors, reqSolvable=reqSolvable, solvable=solvable, reqbackendset=reqbackendset, backendset=backendset, includeMultisession=includeMultisession)
  cph, = executeQuery(conn, query)

  ret = Array{Tuple{Int64,Int64,Symbol},1}()
  @showprogress 1 "Get ExVertex IDs..." for data in cph.results[1]["data"]
    exvid, neoid, vsym = data["row"][1], data["row"][2], Symbol(data["row"][3])
    push!(ret, (exvid,neoid,vsym)  )
  end
  return ret
end

"""
    $(SIGNATURES)

Return array of tuples with ExVertex IDs and Neo4j IDs for vertices with label in session.
"""
function getExVertexNeoIDs(
        conn::Neo4j.Connection;
        label::AS="",
        solvable::Int=1,
        backendset::Int=1,
        session::AS="",
        robot::AS="",
        user::AS="",
        reqbackendset::Bool=true  ) where {AS <: AbstractString}
  #
  sn = length(session) > 0 ? ":"*session : ""
  rn = length(robot) > 0 ? ":"*robot : ""
  un = length(user) > 0 ? ":"*user : ""
  lb = length(label) > 0 ? ":"*label : ""
  query = "match (n$(sn)$(rn)$(un)$(lb)) where n.solvable=$(solvable) and exists(n.exVertexId)"
  query = reqbackendset ? query*" and n.backendset=$(backendset)" : query
  query = query*" return n.exVertexId, id(n), n.label"

  cph, = executeQuery(conn, query)
  # loadtx = transaction(conn)
  # cph = loadtx(query, submit=true)

  ret = Array{Tuple{Int64,Int64,Symbol},1}()
  @showprogress 1 "Get ExVertex IDs..." for data in cph.results[1]["data"]
    exvid, neoid, vsym = data["row"][1], data["row"][2], Symbol(data["row"][3])
    push!(ret, (exvid,neoid,vsym)  )
  end
  return ret
end

"""
    $(SIGNATURES)

Return array of tuples with ExVertex IDs and Neo4j IDs for all poses.
"""
function getPoseExVertexNeoIDs(conn::Neo4j.Connection;
        solvable::Int=1,
        backendset::Int=1,
        session::AS="",
        reqbackendset::Bool=true  ) where {AS <: AbstractString}
  #
  getPoseExVertexNeoIDs(conn,
          label="POSE",
          solvable=solvable,
          backendset=backendset,
          session=session,
          reqbackendset=reqbackendset  )
end

function checkandinsertedges!(fgl::G, exvid::Int, nei::CloudVertex; solvable::Int=1, backendset::Int=1) where G <: AbstractDFG
  if nei.properties["solvable"]==solvable &&
     nei.properties["backendset"] == backendset &&
     haskey(fgl.g.vertices, nei.exVertexId)
    #&& nei.exVertexId <= length(fgl.g.vertices)
    alreadythere = false

    # TODO -- error point
    v2 = fgl.g.vertices[nei.exVertexId]
    for graphsnei in Graphs.out_neighbors(v2, fgl.g) # specifically want Graphs function
      # want to ignore if the edge was previously added from the other side, comparing to the out neighbors in the Graphs structure
      if graphsnei.index == exvid #graphsnei.index == nei.exVertexId
        alreadythere = true
        break;
      end
    end
    if !alreadythere
      # add the edge to graph
      v1 = fgl.g.vertices[exvid]
      makeAddEdge!(fgl, v1, v2, saveedgeID=false)
    end
  end
  nothing
end

function copyAllEdges!(fgl::G, cverts::Dict{Int64, CloudVertex}, IDs::Array{Tuple{Int64,Int64, Symbol},1}) where G <: AbstractDFG
  # TODO -- major opportunity for streamlining and improving performance
  # do entire graph, one node at a time
  @showprogress 1 "Copy all edges..." for ids in IDs
    for nei in CloudGraphs.get_neighbors(fgl.cg, cverts[ids[2]], needdata=true)
      checkandinsertedges!(fgl, ids[1], nei, solvable=1, backendset=1)
    end
  end
  nothing
end

function insertnodefromcv!(fgl::G, cvert::CloudGraphs.CloudVertex) where G <: AbstractDFG
  exvid = cvert.exVertexId
  neoid = cvert.neo4jNodeId
  exvert = cloudVertex2ExVertex(cvert)
  # TODO -- change to addNode
  Graphs.add_vertex!(fgl.g, exvert)
  fgl.id < exvert.index ? fgl.id = exvert.index : nothing
  fgl.cgIDs[exvid] = neoid
  if typeof(exvert.attributes["data"]) == VariableNodeData  # variable node
    fgl.IDs[Symbol(exvert.label)] = exvid
    push!(fgl.nodeIDs, exvid)
  else # function node
    fgl.fIDs[Symbol(exvert.label)] = exvid
    push!(fgl.factorIDs, exvid)
  end
  nothing
end

"""
    $(SIGNATURES)

Copy all variable and factor nodes from DB into fgl as listed in IDs vector.
"""
function copyAllNodes!(fgl::G,
                       cverts::Dict{Int64, CloudVertex},
                       IDs::Array{Tuple{Int64,Int64, Symbol},1}  ) where G <: AbstractDFG
                       # conn::Neo4j.Connection  ) # what the?
  #
  @showprogress 1 "Copy all nodes..." for ids in IDs
    cvert = CloudGraphs.get_vertex(fgl.cg, ids[2], false)
    cverts[ids[2]] = cvert
    insertnodefromcv!(fgl, cvert)
  end
  nothing
end

"""
    $(SIGNATURES)

Copy nodes into `fgl` from DB as listed in the `IDs` vector.
"""
function copyGraphNodesEdges!(fgl::G,
                              IDs::Array{Tuple{Int64,Int64, Symbol},1}  ) where G <: AbstractDFG
  #
  if length(IDs) > 0
    cverts = Dict{Int64, CloudVertex}()
    unsorted = Int64[]
    for ids in IDs
      push!(unsorted, ids[1])
    end
    perm = sortperm(unsorted)

    # get and add all the nodes
    sortedIDs = IDs[perm]
    copyAllNodes!(fgl, cverts, sortedIDs)

    # get and insert all edges
    copyAllEdges!(fgl, cverts, sortedIDs)
    return true
  else
    print(".")
    return false
  end
end

"""
    $(SIGNATURES)

Copy sub graph portion defined by, and including a depth of neighbors, from the lbls vector.
"""
function subLocalGraphCopy!(fgl::G,
                            lbls::Union{Vector{AS}, Vector{Symbol}};
                            neighbors::Int=0,
                            reqbackendset::Bool=true,
                            reqSolvable::Bool=true,
                            includeMultisession::Bool=false) where {G <: AbstractDFG, AS <: AbstractString}
  #
  @warn "subGraphCopy! is a work in progress"
  conn = fgl.cg.neo4j.connection
  IDs = getLblExVertexNeoIDs(conn, string.(lbls), session=fgl.sessionname, robot=fgl.robotname, user=fgl.username, reqbackendset=reqbackendset, reqSolvable=reqSolvable, neighbors=neighbors, includeMultisession=includeMultisession)
  println("fullSubGraphCopy: $(length(IDs)) nodes in session $(fgl.sessionname) if reqbackendset=$reqbackendset and reqSolvable=$reqSolvable...")
  copyGraphNodesEdges!(fgl, IDs)
  nothing
end

"""
    $(SIGNATURES)

Fetch a full copy of the DB factor graph under fgl.sessionname.
"""
function fullLocalGraphCopy!(fgl::G;
                             reqbackendset::Bool=true,
                             reqSolvable::Bool=true  ) where G <: AbstractDFG
  #
  conn = fgl.cg.neo4j.connection
  IDs = getAllExVertexNeoIDs(conn, sessionname=fgl.sessionname, robotname=fgl.robotname, username=fgl.username, reqbackendset=reqbackendset, reqSolvable=reqSolvable)
  println("fullLocalGraphCopy: $(length(IDs)) nodes in subgraph for user=$(fgl.username), robot=$(fgl.robotname), session=$(fgl.sessionname) if reqbackendset=$reqbackendset and reqSolvable=$reqSolvable...")
  copyGraphNodesEdges!(fgl, IDs)
end

"""
    $(SIGNATURES)

Set all Neo4j nodes in this session solvable = 1, warning function does not support new GraffSDK data storage formats.
"""
function setAllDBSolvable!(conn::Neo4j.Connection,
                        sessionname::AS,
                        robotname::AS,
                        username::AS) where {AS <: AbstractString}
  #
  @warn "Obsolete setAllDBSolvable! function, see GraffSDK for example solvable function instead."
  sn = length(sessionname) > 0 ? ":"*sessionname : ""
  rn = length(robotname) > 0 ? ":"*robotname : ""
  un = length(username) > 0 ? ":"*username : ""
  query = "match (n$(sn)$(rn)$(un)) set n.solvable=1"
  cph, loadresult = executeQuery(conn, query)
  nothing
end

# TODO --this will only work with DB version, introduces a bug
function setAllDBSolvable!(fgl::G) where G <: AbstractDFG
  setAllDBSolvable!(fgl.cg.neo4j.connection, fgl.sessionname)
end

function setBackendWorkingSet!(
            conn::Neo4j.Connection,
            sessionname::AbstractString,
            robotname::AbstractString,
            username::AbstractString  )
  #
  sn = length(sessionname) > 0 ? ":"*sessionname : ""
  rn = length(robotname) > 0 ? ":"*robotname : ""
  un = length(username) > 0 ? ":"*username : ""
  query = "match (n$(sn)$(rn)$(un)) where not (n:NEWDATA) set n.backendset=1"

  cph, loadresult = executeQuery(conn, query)
  nothing
end

"""
    $(SIGNATURES)

Obtain Neo4j global database address and login credientials from STDIN, then insert and return in the addrdict colletion.
"""
function askneo4jcredentials!(;addrdict=Dict{AbstractString,AbstractString}() )
  need = ["neo4jHost";"neo4jPort";"neo4jUsername";"neo4jPassword";"session";"robotId"]
  info("Please enter information for Neo4j DB:")
  for n in need
    info(n)
    str = readline(stdin)
    addrdict[n] = str
    if length(str) > 0
      if str[end] == "\n"
        addrdict[n] = str[1:(end-1)]
      end
    end
  end
  return addrdict
end

"""
    $(SIGNATURES)

Obtain Mongo database address and login credientials from STDIN, then insert and return in the addrdict colletion.
"""
function askmongocredentials!(;addrdict=Dict{AbstractString,AbstractString}() )
  need = ["mongoHost";"mongoPort";"mongoUsername";"mongoPassword"]
  info("Please enter information for MongoDB:")
  for n in need
    info(n)
    n == "mongoHost" && haskey(addrdict, "neo4jHost") ? print(string("[",addrdict["neo4jHost"],"]: ")) : nothing
    str = readline(stdin)
    addrdict[n] = str
    if length(str) > 0
      if str[end] == "\n"
        addrdict[n] = str[1:(end-1)]
      end
    end
  end
  addrdict["mongoIsUsingCredentials"] = false,
  if addrdict["mongoHost"] == "" && haskey(addrdict, "neo4jHost")
    addrdict["mongoHost"] = addrdict["neo4jHost"]
  elseif addrdict["mongoHost"] != ""
    nothing
  else
    error("Don't know how to get to MongoDB address.")
  end
  return addrdict
end

"""
    $(SIGNATURES)

Obtain database addresses and login credientials from STDIN, as well as a few case dependent options.
"""
function consoleaskuserfordb(;nparticles=false, drawdepth=false, clearslamindb=false, multisession=false, drawedges=false)
  addrdict = Dict{AbstractString, Any}()  #Union{AbstractString, Vector{String}}
  askneo4jcredentials!(addrdict=addrdict)
  askmongocredentials!(addrdict=addrdict)
  need = String[]
  !nparticles ? nothing : push!(need, "num particles")
  !drawdepth ? nothing : push!(need, "draw depth")
  !clearslamindb ? nothing : push!(need, "clearslamindb")
  !user ? nothing : push!(need, "user")
  !robot ? nothing : push!(need, "robot")
  !multisession ? nothing : push!(need, "multisession")
  !drawedges ? nothing : push!(need, "draw edges")

  info("Please also enter information for:")
  for n in need
    info(n)
    n == "draw depth" ? print("[y]/n: ") : nothing
    n == "draw edges" ? print("[y]/n: ") : nothing
    n == "num particles" ? print("[100]: ") : nothing
    n == "clearslamindb" ? print("yes/[no]: ") : nothing
    n == "user" ? print("[]: ") : nothing
    n == "robot" ? print("[]: ") : nothing
    n == "multisession" ? print("comma separated list session names/[n]: ") : nothing
    str = readline(stdin)
    addrdict[n] = str
    if length(str) > 0
      if str[end] == "\n"
        addrdict[n] = str[1:(end-1)]
      end
    end
  end
  if drawdepth
    addrdict["draw depth"] = addrdict["draw depth"]=="" || addrdict["draw depth"]=="y" || addrdict["draw depth"]=="yes" ? "y" : "n"
  end
  if drawedges
    addrdict["draw edges"] = addrdict["draw edges"]=="" || addrdict["draw edges"]=="y" || addrdict["draw edges"]=="yes" ? "y" : "n"
  end
  if nparticles
    addrdict["num particles"] = addrdict["num particles"]!="" ? addrdict["num particles"] : "100"
  end
  if clearslamindb
    addrdict["clearslamindb"] = addrdict["clearslamindb"]=="" || addrdict["clearslamindb"]=="n" || addrdict["clearslamindb"]=="no" ? "n" : addrdict["clearslamindb"]
  end
  if user
      addrdict["user"] = addrdict["user"]!="" ? addrdict["user"] : ""
  end
  if robot
    addrdict["robot"] = addrdict["robot"]!="" ? addrdict["robot"] : ""
  end
  if multisession
    addrdict["multisession"] = strip.(Vector{String}(split(addrdict["multisession"],',')))
  end
  return addrdict
end

"""
    $(SIGNATURES)

Connect to databases via network according to addrdict, or ask user for credentials and return
active cloudGraph object, as well as addrdict.
"""
function standardcloudgraphsetup(;addrdict=nothing,
            nparticles::Bool=false,
            drawdepth::Bool=false,
            drawedges::Bool=false,
            clearslamindb::Bool=false,
            multisession::Bool=false  )
  #
  if addrdict == nothing
    addrdict = consoleaskuserfordb(nparticles=nparticles, drawdepth=drawdepth, clearslamindb=clearslamindb, multisession=multisession, drawedges=drawedges)
  end

  @warn "Not considering field: addrdict[\"mongoIsUsingCredentials\"]"
  # Connect to database
  # addrdict["mongoIsUsingCredentials"]
  configuration = CloudGraphs.CloudGraphConfiguration(
                             addrdict["neo4jHost"], parse(Int, addrdict["neo4jPort"]), addrdict["neo4jUsername"], addrdict["neo4jPassword"],
                             addrdict["mongoHost"], parse(Int, addrdict["mongoPort"]), false, addrdict["mongoUsername"], addrdict["mongoPassword"]);
  cloudGraph = connect(configuration, IncrementalInference.encodePackedType, Caesar.getpackedtype, IncrementalInference.decodePackedType);
  # conn = cloudGraph.neo4j.connection
  # register types of interest in CloudGraphs
  # registerGeneralVariableTypes!(cloudGraph) # TESTING -- using dispatch instead
  Caesar.usecloudgraphsdatalayer!()

  return cloudGraph, addrdict
end

"""
    $(SIGNATURES)

Walk through vertex bigDataElements and return the last matching description.
"""
function getBigDataElement(vertex::CloudVertex, description::AbstractString)
  bde = nothing
  for bDE in vertex.bigData.dataElements
    if bDE.description == description
      bde = bDE
    end
  end
  return bde
end

"""
    $(SIGNATURES)

Return true if vertex has bigDataElements with matching description.
"""
function hasBigDataElement(vertex::CloudVertex, description::AbstractString)
  for bDE in vertex.bigData.dataElements
    if bDE.description == description
      return true
    end
  end
  return false
end

"""
    $(SIGNATURES)

Append big data element into current blob store and update associated global
vertex information.
"""
function appendvertbigdata!(cloudGraph::CloudGraph,
        cv::CloudVertex,
        description,
        data::Vector{UInt8}  )
  #
  bd = CloudGraphs.read_BigData!(cloudGraph, cv)
  bdei = CloudGraphs.BigDataElement(description, data)
  push!(cv.bigData.dataElements, bdei);
  CloudGraphs.save_BigData!(cloudGraph, cv)
end

"""
    $(SIGNATURES)

Append big data element into current blob store and update associated global
vertex information.
"""
function appendvertbigdata!(fgl::G,
      vert::Graphs.ExVertex,
      description::AbstractString,
      data::Vector{UInt8}  ) where G <: AbstractDFG
  #
  # TODO -- improve get/fetch vertex abstraction
  cvid = fgl.cgIDs[vert.index]
  cv = CloudGraphs.get_vertex(fgl.cg, cvid, true)
  appendvertbigdata!(fgl.cg, cv, description, data)
end

"""
    appendvertbigdata!(fg, sym, descr, data)

Append big data element into current blob store using parent appendvertbigdata!,
but here specified by symbol of variable node in the FactorGraph. Note the
default data layer api definition. User must define dlapi to refetching the
 vertex from the data layer. localapi avoids repeated network database fetches.
"""
function appendvertbigdata!(fgl::G,
      sym::Symbol,
      description::AbstractString,
      data  ) where G <: AbstractDFG
  #
  appendvertbigdata!(fgl,
        getVert(fgl, sym),
        description,
        data  )
end

"""
    $(SIGNATURES)

Fetch and insert list of CloudVertices into FactorGraph object, up to neighbor depth.
"""
function fetchsubgraph!(fgl::G,
                        cvs::Vector{CloudGraphs.CloudVertex};
                        numneighbors::Int=0 ) where G <: AbstractDFG
          # overwrite::Bool=false  )
  # recursion termination condition
  numneighbors >= 0 ? nothing : (return nothing)

  for cv in cvs
    # test if these are already in fgl
    if !hasval(fgl.cgIDs, cv.neo4jNodeId) # may have been inserted as previous neighbor
      # add this vert to graph
      insertnodefromcv!(fgl, cv)

      # recursive call on neighbors here
      neicvs = CloudGraphs.get_neighbors(fgl.cg, cv, needdata=true)
      fetchsubgraph!(fgl, neicvs; numneighbors=numneighbors-1 )
      # add edges associated with the neighbors

      if numneighbors-1 >= 0
        for cvn in neicvs
          checkandinsertedges!(fgl, cv.exVertexId, cvn, solvable=1, backendset=1)
          # makeAddEdge!(fgl, fgl.g.vertices[cv.exVertexId], fgl.g.vertices[cvn.exVertexId], saveedgeID=false)
        end
      end
    end

  end
end

"""
    $SIGNATURES

Fetch and insert list of Neo4j IDs into FactorGraph object, up to neighbor depth.
"""
function fetchsubgraph!(fgl::G,
                        neoids::Vector{Int};
                        numneighbors::Int=0 ) where G <: AbstractDFG
                        # overwrite::Bool=false  )
  #
  for nid in neoids
    # test if these are already in fgl
    if !hasval(fgl.cgIDs, nid)
      cv = CloudGraphs.get_vertex(fgl.cg, nid, false)
      fetchsubgraph!(fgl, [cv], numneighbors=numneighbors )
    end
  end
  nothing
end

"""
    getVertNeoIDs!(::CloudGraph, res::Dict{Symbol, Int}; session::AbstractString="NA", robot::AbstractString="NA", user::AbstractString="NA")

Insert into and return dict `res` with Neo4j IDs of ExVertex labels as stored per session in Neo4j database.
"""
function getVertNeoIDs!(cloudGraph::CloudGraph, res::Dict{Symbol, Int}; session::AbstractString="NA", robot::AbstractString="NA", user::AbstractString="NA")
  loadtx = transaction(cloudGraph.neo4j.connection)
  syms = collect(keys(res))
  query = "match (n:$(session):$robot:$user) where "
  for i in 1:length(syms)
    sym = syms[i]
    query =query*"n.label='$(sym)' "
    if i < length(syms)
      query =query*"or "
    end
  end
  query = query*"return id(n), n.label"
  cph = loadtx(query, submit=true)
  for qr in cph.results[1]["data"]
    res[Symbol(qr["row"][2])] = qr["row"][1]
  end
  # Symbol(cph.results[1]["data"][1]["row"][2]) == sym ? nothing : error("Neo4j query not returning $(sym)")
  return res
end

"""
    removeNeo4jID(cg::CloudGraph, neoid=-1)

Remove node from Neo4j according to Neo4j Node ID. Big data elements that may be associated with this
node are not removed.
"""
function removeNeo4jID(cg::CloudGraph; neoid::Int=-1)
  neoid > 0 ? nothing : error("Can't delete negative neoid=$(neoid).")
  loadtx = transaction(cg.neo4j.connection)
  query =  "match (n) where id(n)=$(neoid) detach delete n return count(n)"
  cph = loadtx(query, submit=true)
  commit(loadtx)
  cnt = Int(cph.results[1]["data"][1]["row"][1])
  cnt == 1 ? nothing : error("Did not delete just one entry after running query = $(query)")
  nothing
end

"""
    getfirstpose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)

Return Tuple{Symbol, Int} of first pose symbol and Neo4j node ID.
"""
function getfirstpose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)
  query = "match (n:$(session):$(robot):$(user):POSE) with n.label as nlbl, n.exVertexId as exvid, id(n) as neoid order by exvid asc limit 1 return nlbl, neoid"
  cph, = executeQuery(cg, query)
  # loadtx = transaction(cg.neo4j.connection)
  # cph = loadtx(query, submit=true)
  Symbol(cph.results[1]["data"][1]["row"][1]), cph.results[1]["data"][1]["row"][2]
end
getFirstPose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString) = getfirstpose(cg, session, robot, user)

"""
    getLastPose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)

Return Tuple{Symbol, Int} of first pose symbol and Neo4j node ID.
"""
function getLastPose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)
  query = "match (n:$(session):$robot:$user:POSE) with n.label as nlbl, n.exVertexId as exvid, id(n) as neoid order by exvid desc limit 1 return nlbl, neoid"
  cph, = executeQuery(cg, query)
  # loadtx = transaction(cg.neo4j.connection)
  # cph = loadtx(query, submit=true)
  Symbol(cph.results[1]["data"][1]["row"][1]), cph.results[1]["data"][1]["row"][2]
end

"""
    insertrobotdatafirstpose!(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString, robotdict::Dict)

Saves robotdict via JSON to first pose in a SESSION in the database. Used for
storing general robot specific data in easily accessible manner. Can fetch later
retrieve same dict with counterpart `fetchrobotdatafirstpose` function.
"""
function insertrobotdatafirstpose!(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString, robotdict::Dict)
  vsym, neoid = getfirstpose(cg, session, robot, user)
  cv = CloudGraphs.get_vertex(cg, neoid, true)
  appendvertbigdata!(cg, cv, "robot_description", json(robotdict).data  )
end

function tryunpackalltypes!(resp::Dict)
  for (k,v) in resp
    tv = typeof(v)
    if tv == Vector{Any}
      ttv = typeof(v[1])
      if ttv == Vector{Any}
        try resp[k] = hcat(v...) catch end
      else
        try resp[k] = Vector{ttv}(v) catch end
      end
    end
  end
  nothing
end

"""
    fetchrobotdatafirstpose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)

Return dict of JSON parsed "robot_description" field as was inserted by counterpart
`insertrobotdatafirstpose!` function. Used for storing general robot specific data
in easily accessible manner.
"""
function fetchrobotdatafirstpose(cg::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)
  vsym, neoid = getfirstpose(cg, session, robot, user)
  cv = CloudGraphs.get_vertex(cg, neoid, true)
  bde = Caesar.getBigDataElement(cv, "robot_description")
  resp = JSON.parse(String(take!(IOBuffer(bde.data))))
  tryunpackalltypes!(resp)
  return resp
end

"""
    getRangeKDEMax2D(cgl::CloudGraph, session::AbstractString, vsym1::Symbol, vsym2::Symbol)

Calculate the cartesian distange between two vertices in the graph, by session and symbol names,
and by maximum belief point.
"""
function getRangeKDEMax2D(cgl::CloudGraph, session::AbstractString, vsym1::Symbol, vsym2::Symbol)
  # get the relavent neo4j ids
  syms = Dict(vsym1=>0, vsym2 => 0)
  getVertNeoIDs!(cgl, syms, session=session)

  # build a local subgraph
  sfg = initfg(cloudgraph=cgl, sessionname=session)
  fetchsubgraph!(sfg, collect(values(syms)), numneighbors=1 )

  # calculate distances on local subgraph
  getRangeKDEMax2D(sfg, vsym1, vsym2)
end

"""
    db2jld(cgl::CloudGraph, session::AbstractString, filename::AbstractString)

Fetch and save a FactorGraph session to a jld, using CloudGraph object and session definition.
"""
function db2jld(cgl::CloudGraph, session::AbstractString, filename::AbstractString)
  fg = initfg(sessionname=session, cloudgraph=cgl)
  fullLocalGraphCopy!(fg)
  savejld(fg, file=filename)
  return fg
end

"""
    db2jld(filename::AbstractString; addrdict::NothingUnion{Dict{AbstractString, AbstractString}}=nothing )

Fetch and save a FactorGraph session to a jld, using or asking STDIN for credentials in the addrdict field.
"""
function db2jld(filename::AbstractString; addrdict::NothingUnion{Dict{AbstractString, AbstractString}}=nothing )
  cg, cr = standardcloudgraphsetup(addrdict=addrdict)
  db2jld(cg, cr["session"], filename)
end

function deleteServerSession!(cloudGraph::CloudGraph, session::AbstractString, robot::AbstractString, user::AbstractString)
  query =  "match (n:$(session):$robot:$user)
            detach delete n
            return count(n)"
  return executeQuery(cloudGraph.neo4j.connection, query)
end

# function syncmongos()
#
# end

  #


text:
# Fusion trees:
#==============================================================================#
struct FusionTree{G<:Sector,N₁,N₂,M,L,T}
    outgoing::NTuple{N₁,G}
    incoming::NTuple{N₂,G}
    innerlines::NTuple{M,G} # M = N₁+N₂-3
    vertices::NTuple{L,T} # L = N₁+N₂-2
end
FusionTree{G,N₁,N₂}(outgoing::NTuple{N₁}, incoming::NTuple{N₂}, innerlines::NTuple{M}, vertices::NTuple{L}) where {G<:Sector,N₁,N₂,M,L} =
    FusionTree{G,N₁,N₂,M,L,vertex_labeltype(G)}(outgoing, incoming, innerlines, vertices)
FusionTree{G,N₁,N₂}(outgoing, incoming, innerlines) where {G<:Sector,N₁,N₂} =
    FusionTree{G,N₁,N₂}(outgoing, incoming, innerlines, ntuple(n->nothing,valsub(valadd(Val(N₁),Val(N₂)),Val(2))))

FusionTree{G}(outgoing::NTuple{N₁}, incoming::NTuple{N₂}, innerlines::NTuple{M}, vertices::NTuple{L}) where {G<:Sector,N₁,N₂,M,L} =
    FusionTree{G,N₁,N₂,M,L,vertex_labeltype(G)}(outgoing, incoming, innerlines, vertices)
FusionTree{G}(outgoing::NTuple{N₁}, incoming::NTuple{N₂}, innerlines) where {G<:Sector,N₁,N₂} =
    FusionTree{G,N₁,N₂}(outgoing, incoming, innerlines, ntuple(n->nothing,valsub(valadd(Val(N₁),Val(N₂)),Val(2))))

FusionTree(outgoing::NTuple{N₁,G}, incoming::NTuple{N₂,G}, innerlines::NTuple{M,G}) where {G<:Sector, N₁, N₂, M} = FusionTree{G,N₁,N₂}(outgoing,incoming,innerlines,ntuple(n->nothing,valsub(valadd(Val(N₁),Val(N₂)),Val(2))))
FusionTree(outgoing::NTuple{N₁,G}, incoming::NTuple{N₂,G}, innerlines::Tuple{}, vertices::Tuple{}) where {G<:Sector, N₁, N₂} = FusionTree{G,N₁,N₂}(outgoing,incoming,innerlines,vertices)

# TODO: don't require innerlines for fusiontype(G) <: Abelian, remove dependence on fusiontrees in _permute(,,,::Type{Abelian})

# Properties
sectortype(::Type{<:FusionTree{G}}) where {G<:Sector} = G
fusiontype(::Type{<:FusionTree{G}}) where {G<:Sector} = fusiontype(G)
numout(::Type{<:FusionTree{G,N₁,N₂}}) where {G<:Sector,N₁,N₂} = N₁
numin(::Type{<:FusionTree{G,N₁,N₂}}) where {G<:Sector,N₁,N₂} = N₂

sectortype(t::FusionTree) = sectortype(typeof(t))
fusiontype(t::FusionTree) = fusiontype(typeof(t))
numout(t::FusionTree) = numout(typeof(t))
numin(t::FusionTree) = numin(typeof(t))

# Fusion tree methods
# adjoint swaps in and out
Base.adjoint(t::FusionTree) = FusionTree(t.incoming, t.outgoing, reverse(t.innerlines), reverse(t.vertices))

Base.split(t::FusionTree{G,0,0}) where {G<:Sector} = (FusionTree{G,0,1}((),(one(G),),()), FusionTree{G,1,0}((one(G),),(),()))
Base.split(t::FusionTree{G,1,0}) where {G<:Sector} = (FusionTree{G,1,1}(t.outgoing,t.outgoing,()), t)
Base.split(t::FusionTree{G,0,1}) where {G<:Sector} = (t, FusionTree{G,1,1}(t.incoming,t.incoming,()))
Base.split(t::FusionTree{G,1,1}) where {G<:Sector} = (t, t)
Base.split(t::FusionTree{G,2,1}) where {G<:Sector} = (t, FusionTree{G,1,1}(t.incoming,t.incoming,()))
Base.split(t::FusionTree{G,1,2}) where {G<:Sector} = (FusionTree{G,1,1}(t.outgoing,t.outgoing,()), t)
function Base.split(t::FusionTree{<:Sector, N₁, N₂}) where {N₁,N₂} # N₁ + N₂ > 3 and N₁,N₂ >= 2
    inner = t.innerlines
    c = inner[N₁-1]
    inner1 = ntuple(n->inner[n], valsub(Val(N₁),Val(2)))
    inner2 = ntuple(n->inner[N₁-1+n], valsub(Val(N₂),Val(2)))
    vert1 = ntuple(n->t.vertices[n], valsub(Val(N₁),Val(1)))
    vert2 = ntuple(n->t.vertices[N₁-1+n], valsub(Val(N₂),Val(1)))
    return FusionTree(t.outgoing, (c,), inner1, vert1), FusionTree((c,), t.incoming, inner2, vert2)
end

function Base.join(t1::FusionTree{G,N₁,N}, t2::FusionTree{G,N,N₂}) where {G<:Sector,N₁,N,N₂}
    t1a, t1b = split(t1)
    t2a, t2b = split(t2)
    t1b == adjoint(t2a) || throw(SectorMismatch())
    t = FusionTree{G,N₁,N₂}(t1a.outgoing, t2b.incoming, (t1a.innerlines..., t1a.incoming[1], t2b.innerlines...), (t1a.vertices..., t2b.vertices...))
    return t=>sqrt(prod(map(dim,t1b.incoming))/dim(t1b.outgoing[1]))
end

# repartition fusion tree
repartition(t::FusionTree, ::Val{N}) where {N} = _repartition(t, Val(N), fusiontype(t))
function _repartition(t::FusionTree{G,N₁,N₂}, ::Val{N}, F::Union{Type{Abelian},Type{SimpleNonAbelian}}) where {G<:Sector,N₁,N₂,N}
    coeff = Bsymbol(one(G), one(G), one(G))
    outer = (t.outgoing..., map(dual, reverse(t.incoming))...)
    innerext = (one(G), first(outer), t.innerlines..., dual(last(outer)), one(G))
    for n = N₁+1:N
         # map fusion vertex c<-(a,b) to splitting vertex (c,dual(b))<-a
        b = dual(outer[n])
        a = innerext[n+1]
        c = innerext[n]
        coeff *= conj(Bsymbol(c, dual(b), a)) # should be a pure phase
    end
    for n = N₁:-1:N+1
        # map splitting vertex (a,b)<-c to fusion vertex a<-(c,dual(b))
        b = outer[n]
        a = innerext[n]
        c = innerext[n+1]
        coeff *= Bsymbol(a,b,c)
    end
    outgoing = ntuple(n->outer[n], Val(N))
    incoming = ntuple(n->dual(outer[N₁+N₂+1-n]), valsub(valadd(Val(N₁),Val(N₂)),Val(N)))
    return ImmutableDict(FusionTree(outgoing,incoming,t.innerlines,t.vertices)=>coeff)
end
# TODO: implementation for DegenerateNonAbelian

# permute fusion tree
function Base.permute(t::FusionTree, p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int}) where {N₁,N₂}
    @assert numin(t) + numout(t) == N₁ + N₂
    @assert braidingtype(sectortype(t)) <: SymmetricBraiding
    return _permute(t::FusionTree, p1, p2, fusiontype(t))
end

function _linearizepermutation(p1::NTuple{N₁,Int}, p2::NTuple{N₂}, n₁::Int, n₂::Int) where {N₁,N₂}
    p1′ = ntuple(Val(N₁)) do n
        p1[n] > n₁ ? n₂+2n₁+1-p1[n] : p1[n]
    end
    p2′ = ntuple(Val(N₂)) do n
        p2[N₂+1-n] > n₁ ? n₂+2n₁+1-p2[N₂+1-n] : p2[N₂+1-n]
    end
    return (p1′..., p2′...)
end

function _permute(tsrc::FusionTree, p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int}, ::Type{Abelian}) where {N₁,N₂}
    p = _linearizepermutation(p1, p2, numout(tsrc), numin(tsrc))
    t, coeff = first(repartition(tsrc, valadd(Val(N₁),Val(N₂))))
    for i = 1:(N₁+N₂)
        for j = 1:i-1
            if p[j] > p[i]
                a, b = t.outgoing[j], t.outgoing[i]
                c = first(a ⊗ b)
                coeff *= Rsymbol(a, b, c)
            end
        end
    end
    t = first(fusiontrees(tpermute(t.outgoing, p),()))
    tdst, coeff2 = first(repartition(t, Val(N₁)))
    return ImmutableDict(tdst=>(coeff*coeff2))
end

function _permute(tsrc::FusionTree, p1::NTuple{N₁,Int}, p2::NTuple{N₂,Int}, ::Type{<:NonAbelian}) where {N₁,N₂}
    p = _linearizepermutation(p1, p2, numout(tsrc), numin(tsrc))
    trees = Dict(t=>c for (t,c) in repartition(tsrc, valadd(Val(N₁),Val(N₂))))
    newtrees = similar(trees)
    for s in permutation2swaps(p)
        for (t, c) in trees
            for (t′,c′) in swap(t, s)
                newtrees[t′] = get(newtrees, t′, zero(c)) + c*c′
            end
        end
        trees, newtrees = newtrees, trees
        empty!(newtrees)
    end
    s = start(trees)
    (t, c), s = next(trees, s)
    tlist = repartition(t, Val(N₁))
    output = Dict{keytype(tlist),valtype(tlist)}()
    for (t′, c′) in tlist
        output[t′] = c*c′
    end
    while !done(trees, s)
        v, s = next(trees, s)
        t, c = v
        for (t′,c′) in repartition(t, Val(N₁))
            output[t′] = get(output, t′, zero(c)) + c*c′
        end
    end
    return output
end

"""
    function swap(t::FusionTree{<:Sector,N,0}, i) where {N} -> ImmutableDict{typeof(t),<:Number}

Performs a swap of neighbouring outgoing indices `i` and `i+1` on a fusion tree `t`,
and returns the result as a linked list of output trees and corresponding coefficients.
"""
function swap(t::FusionTree{<:Sector,N,0}, i) where {N}
    1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
    _swap(t, i, fusiontype(t))
end

function _swap(t::FusionTree{<:Sector,N,0}, i, ::Type{Abelian}) where {N}
    outer = t.outgoing
    inner = t.innerlines
    b, d = outer[i], outer[i+1]
    outer = setindex(outer, d, i)
    outer = setindex(outer, b, i+1)
    if i > 1 && i < N-1
        e = (i == N-2 ? dual(outer[N]) : inner[i])
        inner = setindex(inner, first(e ⊗ dual(b)), i-1)
    end
    return ImmutableDict(FusionTree(outer, (), inner) => Rsymbol(b, d, first(b ⊗ d)))
end

function _swap(t::FusionTree{<:Sector,N,0}, i, ::Type{SimpleNonAbelian}) where {N}
    outer = t.outgoing
    inner = t.innerlines
    b = outer[i]
    d = outer[i+1]
    a = (i == 1 ? one(b) : (i == 2 ? outer[1] : inner[i-2]))
    c = (i == 1 ? b : (i == N-1 ? dual(d) : inner[i-1]))
    e = (i == N-1 ? one(b) : (i == N-2 ? dual(outer[N]) : inner[i]))
    outer = setindex(outer, d, i)
    outer = setindex(outer, b, i+1)
    iter = a ⊗ d
    s = start(iter)
    c′, s = next(iter, s)
    while !Nsymbol(b, c′, e)
        c′, s = next(iter, s)
    end
    coeff = conj(Rsymbol(b,a,c))*Fsymbol(b,a,d,e,c,c′)*Rsymbol(b,c′,e)
    inner′ = 1 < i < N-1 ? setindex(inner, c′, i-1) : inner
    output = ImmutableDict(FusionTree(outer, t.incoming, inner′, t.vertices)=>coeff)
    while !done(iter, s)
        c′, s = next(iter, s)
        Nsymbol(b, c′, e) || continue
        coeff = conj(Rsymbol(b,a,c))*Fsymbol(b,a,d,e,c,c′)*Rsymbol(b,c′,e)
        inner′ = 1 < i < N-1 ? setindex(inner, c′, i-1) : inner
        output = ImmutableDict(output, FusionTree(outer, t.incoming, inner′, t.vertices)=>coeff)
    end
    return output
end

# TODO: implementation for DegenerateNonAbelian

# Fusion tree iterators
include("iterator.jl")

# Show methods
function Base.show(io::IO, t::FusionTree{G,0,0}) where {G<:Sector}
    if get(io, :compact, false)
        println(io, "FusionTree{", G, ", 0, 0}((), (), (), ())")
    else
        println(io, "FusionTree{", G, ", 0, 0}: (empty fusion tree)")
    end
end
function Base.show(io::IO, t::FusionTree{G,1,0}) where {G<:Sector}
    if get(io, :compact, false)
        println(io, "FusionTree{", G, ", 1, 0}(", t.outgoing ,"(), (), ())")
    else
        println(io, "FusionTree{", G, ", 1, 0}: (vacuum line)")
    end
end
function Base.show(io::IO, t::FusionTree{G,0,1}) where {G<:Sector}
    if get(io, :compact, false)
        println(io, "FusionTree{", G, ", 0, 1}((), ", t.incoming ,"(), ())")
    else
        println(io, "FusionTree{", G, ", 0, 1}: (dual vacuum line)")
    end
end

function Base.show(io::IO, t::FusionTree{G,N₁,N₂,M,K,Void}) where {G<:Sector,N₁,N₂,M,K}
    if get(io, :compact, false)
        print(io, "FusionTree{", G, ", ", N₁, ", ", N₂, "}(", t.outgoing, ", ", t.incoming, ", ", t.innerlines, ")")
        return
    end
    println(io, "FusionTree{", G, ", ", N₁, ", ", N₂, "}:")
    cprint = x->sprint(showcompact, x)
    outlabels = map(cprint, t.outgoing)
    inlabels = map(cprint, t.incoming)
    L = max(map(strwidth, outlabels)..., map(strwidth, inlabels)...)
    up = lpad(rpad(" ∧ ", L>>1), L)
    down = lpad(rpad(" ∨ ", L>>1), L)
    L = strwidth(up)
    f = lpad(rpad("⋅+-",L>>1,"-"), L, " ")
    l = lpad(rpad("-+⋅",L>>1," "), L, "-")
    m = lpad(rpad("-+-",L>>1,"-"), L, "-")
    a = "-<-"
    pad = " "^L
    line1 = ""
    line2 = ""
    line3 = ""
    line4 = ""
    n = 1
    while n <= N₁ + N₂
        # print incoming or outgoing label
        if n <= N₁
            label = lpad(rpad(outlabels[n], 2), L)
            line4 *= label
            line3 *= down
        else
            label = lpad(rpad(inlabels[N₁+N₂+1-n], 2), L)
            line4 *= label
            line3 *= up
        end
        if n == 1
            line2 *= f
        elseif n == N₁+N₂
            line2 *= l
        else
            line2 *= m
        end
        line1 *= pad
        # print inner label
        if n == 1
            label = N₁ > 0 ? cprint(t.outgoing[1]) : cprint(dual(t.incoming[N₂]))
            label = lpad(rpad(label, 2), strwidth(a))
            line1 *= label
            ℓ = strwidth(label)
            line2 *= lpad(rpad(a, ℓ>>1, "-"), ℓ, "-")
            line3 *= " "^ℓ
            line4 *= " "^ℓ
        elseif n == N₁+N₂-1
            label = N₂ > 0 ? cprint(t.incoming[1]) : cprint(dual(t.outgoing[N₁]))
            label = lpad(rpad(label, 2), strwidth(a))
            line1 *= label
            ℓ = strwidth(label)
            line2 *= lpad(rpad(a, ℓ>>1, "-"), ℓ, "-")
            line3 *= " "^ℓ
            line4 *= " "^ℓ
        elseif n != N₁+N₂
            label = cprint(t.innerlines[n-1])
            label = lpad(rpad(label, 2), strwidth(a))
            line1 *= label
            ℓ = strwidth(label)
            line2 *= lpad(rpad(a, ℓ>>1, "-"), ℓ, "-")
            line3 *= " "^ℓ
            line4 *= " "^ℓ
        end
        n += 1
    end
    println(io, line1)
    println(io, line2)
    println(io, line3)
    println(io, line4)
end

# auxiliary routines

# _abelianinner: generate the inner indices for given outer indices in the abelian case
_abelianinner(outer::Tuple{}) = ()
_abelianinner(outer::Tuple{G}) where {G<:Sector} = ()
_abelianinner(outer::Tuple{G,G}) where {G<:Sector} = ()
_abelianinner(outer::Tuple{G,G,G}) where {G<:Sector} = ()
function _abelianinner(outer::NTuple{N,G}) where {G<:Sector,N}
    c = first(outer[1] ⊗ outer[2])
    return (c, _abelianinner((c, tail2(outer)...)))
end

ext: jl


text:
"
    MPSComoving(leftstate,window,rightstate)

    muteable window of tensors on top of an infinite chain
"
mutable struct MPSComoving{A<:GenericMPSTensor,B<:MPSBondTensor} <: AbstractMPS
    left_gs::InfiniteMPS{A,B}

    ALs::Vector{Union{Missing,A}}
    ARs::Vector{Union{Missing,A}}
    ACs::Vector{Union{Missing,A}}
    CLs::Vector{Union{Missing,B}}

    right_gs::InfiniteMPS{A,B}

    function MPSComoving{A,B}(left,ALs::Vector{Union{Missing,A}},
                            ARs::Vector{Union{Missing,A}},
                            ACs::Vector{Union{Missing,A}},
                            CLs::Vector{Union{Missing,B}},right) where {A<:GenericMPSTensor,B<:MPSBondTensor}
        #todo:insert checks
        return new{A,B}(left,ALs,ARs,ACs,CLs,right);
    end
end

#allow construction with an array of tensors
function MPSComoving(left::InfiniteMPS{A,B},site_tensors::Array{A},right::InfiniteMPS{A,B}) where {A<:GenericMPSTensor,B<:MPSBondTensor}
    site_tensors = copy.(site_tensors)
    for i in 1:length(site_tensors)-1
        (site_tensors[i],C) = leftorth!(site_tensors[i],alg=QRpos());
        site_tensors[i+1] = _permute_front(C*_permute_tail(site_tensors[i+1]))
    end

    (site_tensors[end],C) = leftorth!(site_tensors[end],alg=QRpos());

    CLs = Vector{Union{Missing,B}}(missing,length(site_tensors)+1)
    ALs = Vector{Union{Missing,A}}(missing,length(site_tensors))
    ARs = Vector{Union{Missing,A}}(missing,length(site_tensors))
    ACs = Vector{Union{Missing,A}}(missing,length(site_tensors))

    ALs.= site_tensors;
    CLs[end] = C;

    MPSComoving{A,B}(left,ALs,ARs,ACs,CLs,right)
end
#todo : add constructor given an infinitemps + length

# allow construction with one large tensorkit space
function MPSComoving(f, elt,P::ProductSpace, args...; kwargs...)
    return MPSComoving(f, elt, collect(P), args...; kwargs...)
end

# allow construction given only a physical space and length
function MPSComoving(f,elt, N::Int, V::VectorSpace, args...; kwargs...)
    return MPSComoving(f, elt,fill(V, N), args...; kwargs...)
end

function MPSComoving(f, elt, physspaces::Vector{<:Union{S,CompositeSpace{S}}}, maxvirtspace::S,
                    leftgs::M,
                    rightgs::M) where {S<:ElementarySpace,M<:InfiniteMPS}

    left = virtualspace(leftgs,0);
    right = virtualspace(rightgs,length(physspaces));

    N = length(physspaces)
    virtspaces = Vector{S}(undef, N+1)
    virtspaces[1] = left
    for k = 2:N
        virtspaces[k] = infimum(fuse(virtspaces[k-1], fuse(physspaces[k])), maxvirtspace)
    end
    virtspaces[N+1] = right
    for k = N:-1:2
        virtspaces[k] = infimum(virtspaces[k], fuse(virtspaces[k+1], flip(fuse(physspaces[k]))))
    end
    return MPSComoving(f, elt,physspaces, virtspaces,leftgs,rightgs)
end
function MPSComoving(f,elt,
                    physspaces::Vector{<:Union{S,CompositeSpace{S}}},
                    virtspaces::Vector{S},
                    leftgs::M,
                    rightgs::M) where {S<:ElementarySpace,M<:InfiniteMPS}
    N = length(physspaces)
    length(virtspaces) == N+1 || throw(DimensionMismatch())

    tensors = [TensorMap(f, elt,virtspaces[n] ⊗ physspaces[n], virtspaces[n+1]) for n=1:N]

    return MPSComoving(leftgs,tensors,rightgs)
end

#take a window from an infinite mps, and use that as state
function MPSComoving(state::InfiniteMPS{A,B},len::Int) where {A,B}
    CLs = Vector{Union{Missing,B}}(missing,len+1)
    ALs = Vector{Union{Missing,A}}(missing,len)
    ARs = Vector{Union{Missing,A}}(missing,len)
    ACs = Vector{Union{Missing,A}}(missing,len)

    ALs.= state.AL[1:len];
    ARs.= state.AR[1:len];
    ACs.= state.AC[1:len];
    CLs.= state.CR[0:len];

    MPSComoving{A,B}(state,ALs,ARs,ACs,CLs,state)
end

Base.copy(state::MPSComoving{A,B}) where {A,B} = MPSComoving{A,B}(state.left_gs,copy(state.ALs),copy(state.ARs),copy(state.ACs),copy(state.CLs),state.right_gs);

Base.length(state::MPSComoving) = length(state.ALs)
Base.size(psi::MPSComoving, i...) = size(psi.ALs, i...)

Base.eltype(::Type{MPSComoving{Mtype,Vtype}}) where {Mtype<:GenericMPSTensor,Vtype<:MPSBondTensor} = Mtype
bond_type(::Type{MPSComoving{Mtype,Vtype}}) where {Mtype<:GenericMPSTensor,Vtype<:MPSBondTensor} = Vtype

TensorKit.space(psi::MPSComoving{<:MPSTensor}, n::Integer) = space(psi.AC[n], 2)
virtualspace(psi::MPSComoving, n::Integer) =
    n < length(psi) ? _firstspace(psi.AC[n+1]) : dual(_lastspace(psi.AC[n]))

r_RR(state::MPSComoving)=r_RR(state.right_gs,length(state))
l_LL(state::MPSComoving)=l_LL(state.left_gs,1)

function Base.getproperty(psi::MPSComoving,prop::Symbol)
    if prop == :AL
        return ALView(psi)
    elseif prop == :AR
        return ARView(psi)
    elseif prop == :AC
        return ACView(psi)
    elseif prop == :CR
        return CRView(psi)
    else
        return getfield(psi,prop)
    end
end

function max_Ds(f::MPSComoving{G}) where G<:GenericMPSTensor{S,N} where {S,N}
    Ds = [dim(space(f.left_gs.AL[1],1)) for v in 1:length(f)+1];
    for i in 1:length(f)
        Ds[i+1] = Ds[i]*prod(map(x->dim(space(f.AC[i],x)),ntuple(x->x+1,Val{N-1}())))
    end

    Ds[end] = dim(space(f.right_gs.AL[1],1));
    for i in length(f):-1:1
        Ds[i] = min(Ds[i],Ds[i+1]*prod(map(x->dim(space(f.AC[i],x)),ntuple(x->x+1,Val{N-1}()))))
    end
    Ds
end
function Base.:*(psi::MPSComoving, a::Number)
    return rmul!(copy(psi),a)
end

function Base.:*(a::Number, psi::MPSComoving)
    return lmul!(a,copy(psi))
end

function TensorKit.lmul!(a::Number, psi::MPSComoving)
    psi.ACs .*=a;
    psi.CLs .*=a;
    return psi
end

function TensorKit.rmul!(psi::MPSComoving,a::Number)
    psi.ACs .*=a;
    psi.CLs .*=a;
    return psi
end

function TensorKit.dot(psi1::MPSComoving, psi2::MPSComoving)
    length(psi1) == length(psi2) || throw(ArgumentError("MPS with different length"))
    psi1.left_gs == psi2.left_gs || throw(ArgumentError("left InfiniteMPS is different"))
    psi1.right_gs == psi2.right_gs || throw(ArgumentError("right InfiniteMPS is different"))

    ρr = transfer_right(r_RR(psi2),psi2.AR[2:end],psi1.AR[2:end]);
    return tr(_permute_front(psi1.AC[1])' * _permute_front(psi2.AC[1]) * ρr)
end

TensorKit.norm(psi::MPSComoving) = norm(psi.AC[1])

TensorKit.normalize!(psi::MPSComoving) = rmul!(psi, 1/norm(psi))
TensorKit.normalize(psi::MPSComoving) = normalize!(copy(psi))


text:
module Raylib

export Game, RayColor, WHITE, LIGHTGREY,
  clear_background, draw_text, draw_circle, is_key_down,
  Keys

using Raylib_jll

module Keys
include("Keys.jl")
end

struct RayColor
  r::UInt8
  g::UInt8
  b::UInt8
  a::UInt8
end

const WHITE = RayColor(255,255,255,255)
const LIGHTGREY = RayColor(180,180,180,255)

abstract type Game end

struct WindowSettings
  width::Int
  height::Int
  title::String
  fps::Int
end

width(g::Game) = 800

height(g::Game) = 450

title(g::Game) = "raylib game in Julia"

fps(g::Game) = 60

settings(g::Game) = WindowSettings(width(g), height(g), title(g), fps(g))

draw(g::Game) = nothing

update!(g::Game) = nothing

function Base.run(game::Game)
  s = settings(game)
  ccall((:InitWindow, libraylib), Cvoid, (Cint,Cint,Cstring), s.width, s.height, s.title)
  ccall((:SetTargetFPS, libraylib), Cvoid, (Cint,), s.fps)
  while ccall((:WindowShouldClose, libraylib), Cuchar, ()) == 0
    update!(game)
    ccall((:BeginDrawing, libraylib), Cvoid, ())
    draw(game)
    ccall((:EndDrawing, libraylib), Cvoid, ())
  end
  ccall((:CloseWindow, libraylib), Cvoid, ())
end

clear_background(c::RayColor) =
  ccall((:ClearBackground, libraylib), Cvoid, (RayColor,), c)

draw_text(t::String, x::Int, y::Int, w::Int, c::RayColor) =
  ccall((:DrawText, libraylib), Cvoid, (Cstring, Cint, Cint, Cint, RayColor), t, x, y, w, c)

struct RayVector2
  x::Float32
  y::Float32
  RayVector2(v::AbstractVector{<:Real}) = new(convert(Float32, v[1]),convert(Float32, v[2]))
end

draw_circle(v::AbstractVector{<:Real}, r::Real, c::RayColor) =
  ccall((:DrawCircleV, libraylib), Cvoid, (RayVector2, Cfloat, RayColor),
        RayVector2(v), convert(Float32, r), c)

is_key_down(k::Int) =
  ccall((:IsKeyDown, libraylib), Cuchar, (Cint,), k) == 1

end


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "LOVO21"
cases = [""]
SUITE[model] = BenchmarkGroup()

include("lotka_volterra.jl")
validation = []
final_area = []
intersect_time = []

# ----------------------------------------
#  Case 1:
# ----------------------------------------
const num_divs = 32
const X0minced = unique(mince(X0, num_divs))

# warm-up run
lotka_volterra(X0) = @ivp(x'= lotka_volterra!(x), dim:2, x(0) ∈ X0)
prob = lotka_volterra(X0)
alg = TMJets21a(abstol=1e-10, orderT=7, orderQ=1);
# alg = TMJets21b(abstol=1e-23, orderT=18, orderQ=1)
sol = solve(prob, T=T_lv, alg=alg);

prob_minced = @. lotka_volterra(X0minced);
telapsed = @elapsed sols_minced = @. solve(prob_minced, T=T_lv, alg=alg);
@show(telapsed)

# Results: get times not spent stricktly outside and areas
res = lv_property.(sols_minced, r₀)
time1 = diam(hull(getindex.(res, 1)...))
time2 = diam(hull(getindex.(res, 2)...))
areas = sum(getindex.(res, 3))
ib1 = getindex.(res, 4)
ib2 = getindex.(res, 5)

push!(validation, Int(true))
push!(final_area, trunc(areas, sigdigits=3))
push!(intersect_time, trunc(time1, sigdigits=3), trunc(time2, sigdigits=3))
println("Final area : $(areas)")
println("Time spent in guard : $(time1), $(time2)")

# benchmark
# SUITE[model][cases[1]] = @benchmarkable :(@. solve($prob_minced, T=$T_lv, alg=$alg))
# SUITE[model][cases[1]] = @benchmarkable solve.($prob, T = $T_lv, alg = $alg)

# const ε_ext = 1e-4
# const r₀ = 0.161
# prob = lotka_volterra(; nsplit=7, ε_ext=ε_ext, r₀=r₀)
# alg = TMJets(abstol=1e-14, orderT=7, orderQ=1)

# # warm-up run
# sol_lv = solve(prob, T=T_lv,
#                   alg=alg,
#                   max_jumps=100,
#                   intersect_source_invariant=false,
#                   intersection_method=BoxIntersection(),
#                   clustering_method=BoxClustering(3),
#                   disjointness_method=BoxEnclosure());
# solz_lv = overapproximate(sol_lv, Zonotope);

# # obtain area
# area, time_in_guard = lv_property(solz_lv, ε_ext, r₀)
# push!(validation, Int(true))
# push!(final_area, trunc(area, sigdigits=3))
# push!(intersect_time, trunc(time_in_guard, sigdigits=3))
# println("Final area : $(area)")
# println("Time spent in guard : $(time_in_guard)")

# # benchmark
# SUITE[model][cases[1]] = @benchmarkable solve($prob,
#                   T = $T_lv,
#                   alg = $alg,
#                   max_jumps = 100,
#                   intersect_source_invariant = false,
#                   intersection_method = BoxIntersection(),
#                   clustering_method = BoxClustering(3),
#                   disjointness_method = BoxEnclosure())

# # ==============================================================================
# # Execute benchmarks and save benchmark results
# # ==============================================================================

# # tune parameters
# tune!(SUITE)

# # run the benchmarks
# results = run(SUITE, verbose=true)

# # return the sample with the smallest time value in each test
# println("minimum time for each benchmark:\n", minimum(results))

# # return the median for each test
# println("median time for each benchmark:\n", median(results))

# # export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
    # t = median(results[model][c]).time * 1e-9
    runtimes[c] = telapsed
end

for (i, c) in enumerate(cases)
    print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c])," *
        " $(final_area[i]), $(intersect_time[i])\n")
end

# ==============================================================================
# Create plots
# ==============================================================================
const θ = range(0.0, 2*pi, length=1000)

fig = plot(
    xlims=(0.6, 1.45), ylims=(0.6, 1.4), aspect_ratio=:equal, legend=:none,
    palette = :seaborn_colorblind,
    tickfont=font(30, "Times"), guidefontsize=45,
    xtick=0.6:0.2:1.4, ytick=0.6:0.2:1.4,
    bottom_margin=6mm, left_margin=2mm, right_margin=8mm, top_margin=3mm,
    size=(1000, 1000))

# Flowpipes
for ind in 1:num_divs
    plot!(fig, sols_minced[ind], vars=(1,2), alpha=1.0, color=:blue, linecolor=:match)
end
# Flowpipes crossing guards
for ind in 1:num_divs
    plot!(fig, ib1[ind], alpha=0.9, color=:lightgreen, linecolor=:match)
    plot!(fig, ib2[ind], alpha=0.9, color=:orange, linecolor=:match)
end
plot!(fig, 1 .+ r₀.*cos.(θ), 1 .+ r₀.*sin.(θ), color=:red, label=:none)

xlabel!(fig, L"x")
ylabel!(fig, L"y")

# outside_idx = findall(x -> x == 1, location.(sol_lv))
# inside_idx = findall(x -> x == 2 || x == 3, location.(sol_lv))

# for i in outside_idx
#     plot!(fig, solz_lv[i], vars=(1, 2), lw=0.0, alpha=1.0, color=:blue)
# end

# for i in inside_idx
#     plot!(fig, solz_lv[i], vars=(1, 2), lw=0.0, alpha=1.0, color=:lightgreen)
# end

# B = Ball2([1.0, 1.0], 0.15) # "exact"
# B_ext = overapproximate(B, ε_ext) # outer approximation
# plot!(fig, B, 1e-4, color=:white, lw=2.0, linecolor=:red, tickfont=font(30, "Times"),
#         guidefontsize=45,
#         xlab=L"x",
#         ylab=L"y",
#         xtick=[0.8, 1.0, 1.2], ytick=[0.6, 0.8, 1.0, 1.2],
#         xlims=(0.6, 1.4), ylims=(0.6, 1.4),
#         bottom_margin=6mm, left_margin=2mm, right_margin=8mm, top_margin=3mm,
#         size=(1000, 1000))

savefig("ARCH-COMP21-JuliaReach-LotkaVolterra.png")


text:
include("ThetisSupport.jl")
ThetisSupport.bar(convert(Cint,1))

max_stars_repo_name: AkillesAILimited/Thetis

max_issues_repo_name: AkillesAILimited/Thetis

max_forks_repo_name: AkillesAILimited/Thetis


text:
module Slide

"Slide"
Slide

using Base.Markdown
export pages

pages = [

  md"""
  # module JuliaIntro
  """

  md"""
  # Julia
    $(Dates.today())
  """

    md"""
    # Julia
     - http://julialang.org/
     - 짤 logo_hires.png
    """

  md"""
  # JuliaCon
   - http://juliacon.org/
  """

    md"""
    # JuliaCon 2014
     - http://juliacon.org/2014/
    """

    md"""
    # JuliaCon 2015
     - http://juliacon.org/2015/
    """

  md"""
  # 개발
  """

    md"""
    # JuliaBox
     - https://juliabox.org/
    """

    md"""
    # REPL
    ```
    $ julia
    ```
    """

    md"""
    # IJulia
     - https://github.com/JuliaLang/IJulia.jl
    """

    md"""
    # Juno
     - http://junolab.org/
    """

  md"""
  # 문서
  """

    md"""
    # 책, 튜토리얼, 영상
     - http://julialang.org/learning/
     - 짤 books.png
    """

    md"""
    # 한글 문서
     - https://juliakorea.github.io/latest/
     - 줄리아 언어를 소개합니다 (The Julia Language) BY 권창현 http://thoughts.chkwon.net/the-julia-language/
    """

    md"""
    # 코드
     - https://github.com/JuliaLang/julia

    ```shell
    ~/head/julia master$ git log -p
    ```
    """

  md"""
  # 한국 커뮤니티
  """

    md"""
    # 페이스북
     - https://www.facebook.com/groups/juliakorea
    """

    md"""
    # 슬랙
     - https://juliakorea.slack.com
    """

  md"""
  # 번역 프로세스
  """

    md"""
    # 한글 문서
     - https://juliakorea.github.io/latest/
    """

    md"""
    # doc
     - https://github.com/juliakorea/doc
     - https://github.com/juliakorea/doc/graphs/contributors
    """

    md"""
    # git 도우미
     - https://github.com/juliakorea/doc/wiki/git
    """

    md"""
    # doc 이슈
     - https://github.com/juliakorea/doc/issues
    """

    md"""
    # doc 풀 리퀘스트
     - https://github.com/juliakorea/doc/pulls
    """

    md"""
    # slack 로그
     - https://juliakorea.slack.com/messages/translate/
     - https://juliakorea.slack.com/messages/log/
    """

  md"""
  # 꾸준한 코드 관리, 커버리지
  """

    md"""
    # Travis CI
     - https://travis-ci.org/juliakorea/doc
    """

    md"""
    # AppVeyor
     - https://ci.appveyor.com/project/wookay/doc
    """

    md"""
    # Coveralls
     - https://coveralls.io/github/juliakorea/doc?branch=master
    """

  md"""
  # 실무용 언어
  """

    md"""
    # julia-jobs 메일링
     - https://groups.google.com/forum/#!forum/julia-jobs
    """

    md"""
    # Julia Computing
     - http://juliacomputing.com/
     - http://venturebeat.com/2015/05/18/why-the-creators-of-the-julia-programming-language-just-launched-a-startup/
     - 짤 jcl-logo.png
     - 짤 creators.png
    """

  md"""
  # 연구용 언어
  """

    md"""
    # teaching
     - http://julialang.org/teaching/
    """

    md"""
    # 논문들
     - http://julialang.org/publications/
    """

    md"""
    # Julia Summer of Code 2015
     - http://julialang.org/jsoc
    """

  md"""
  # 과학
  """

    md"""
    # 커뮤니티
     - http://julialang.org/community/

    * JuliaLang – The language
    * JuliaStats – Statistics
    * JuliaOpt – Optimization
    * JuliaParallel – Parallel programming in Julia
    * BioJulia - Biology
    * JuliaAstro – Astronomy
    * JuliaQuant – Finance
    * JuliaQuantum – Julia libraries for quantum science and technology
    * JuliaSparse – Sparse matrix solvers
    * JuliaDiff – Differentiation tools
    * JuliaGPU – GPU computing
    * JuliaWeb – Web stack
    * JuliaCI – Continuous Integration Support for Julia packages
    * JuliaDSP – Digital signal processing
    * JuliaGraphs – Graph Theory and Implementation
    * JuliaLangEs – Julia resources in the Spanish language
    """

  md"""
  # 패키지
   - http://pkg.julialang.org/
  """

    md"""
    # Cxx.jl
     - https://github.com/Keno/Cxx.jl
     - https://github.com/juliakorea/OpenCV.jl/wiki/using-Cxx.jl
    """

    md"""
    ```julia
    using Cxx

    cxx\"\"\"
    // http://todayhumor.com/?programmer_12824
    void navs() {
      int i, j;
      printf("-----\n");
      for (i = 0x00477dc4, j = 0x01084210; i; i >>= 1, j >>= 1) {
        printf("%c%c", (i&1)*('*'-32) + 32, (j&1)*10);
      }
      printf("\n");
      printf("-----\n");
    }
    \"\"\"

    @cxx navs()
    ```
    """

    md"""
    # 유용한 패키지들
     - https://github.com/JuliaOpt/Convex.jl
     - https://github.com/JuliaOpt/JuMP.jl
    """

  md"""
  # 특징.jl
  """

    md"""
    # PARAMETRIC TYPES
    ```julia
    type 만세{짱짱<:Int}
    end

    만세{Int}
    # 만세{String}
    # 만세{Any}

    만세{Int} <: 만세
    # false == (만세{Int} <: Int)
    ```
    """

    md"""
    # PARAMETRIC TYPES - Array
    ```julia
    Array{Int}([1 2 3;4 5 6])
    Array{String}(["밥" "벌이"])
    Array{Any}([1 2 "밥"])
    # Array{Int}(["밥"])
    ```
    """

    md"""
    # PARAMETRIC TYPES - Dict
    ```julia
    Dict{Int,Int}(1=>2, 2=>3)
    Dict{Any,Int}(1=>2, "밥"=>3)
    # Dict{Int,Int}(1=>2, "밥"=>3)
    ```
    """

    md"""
    # PARAMETRIC TYPES - Union
    ```julia
    Union{Int,String}("밥")
    Union{Int,String}(1)
    ```
    """

    md"""
    # MULTIPLE DISPATCH
    ```julia
    멀() = 1
    멀(::Int) = 2
    멀(::Function) = 3

    멀() + 멀(0) + 멀(멀)
    # 멀("")

    methods(멀)
    Base.return_types(멀, (Any,))
    Base.return_types(멀, ())
    ```
    """

    md"""
    # MACROS
    ```julia
    macro 매크로()
      println("안녕")
    end

    @매크로
    ```
    """

    md"""
    # MACROS - PREFIXED STRING LITERALS
    ```julia
    macro 세번_str(a)
      [a,a,a]
    end

    세번"굳"

    v"1.2.3"
    v"1.2.3".major
    v"1.2.3".minor
    v"1.2.3".patch
    ```
    """

    md"""
    # 유용한 개념들
     - type inference
     - type safety
     - abstract, concrete
     - mutable, immutable
     - Any, Base.Bottom
     - linear algebra
       - matrix, vector, scalar
       - range, interval
       - dense, sparse
     - meta programming
     - staged programming
    """
    # - invariant, covariant
    # - parallel tasks, distributed arrays
    # - stream, pipe, async, sync

  md"""
  # 소소한 프로그래머
  """

    md"""
    # 줄리아로 짜본거
     - https://github.com/wookay/journal/wiki/줄리아로짠거
    """

    md"""
    # 줄리아는 뭐다
     - 새로 나온 게임
     - 정원 가꾸기
     - 설계자들 놀이터
     - 청년 혁명
    """

  #= md"""
  # 프로그래밍 스타일
   - 고민하는 지점
  """

    md"""
    # 자유도
     - 표현력
     - 재미
    """

    md"""
    # 팀
     - 줄리아 짱짱맨들
    """

    md"""
    # 패러다임
     - 변화와 수용
     - 생각의 틀
     - 생각의 한계
    """

    md"""
    # 권력과 현실
     - 자본
     - 주인
     - 자주성

       - 만들어야 할 가치
       - 지켜야 할 가치
    """

    md"""
    # 機械時代
     - 학습과 공유
     - 통섭과 환원
     - 성장과 생존
    """
  =#

  md"""
  # 커밍순
  """

    md"""
    # 마일스톤
     - https://github.com/julialang/julia/milestones
    """

    md"""
    # JuliaCon 2016
     - http://juliacon.org/
    """

  md"""
  # 사진 출처
   - Julia
     - logo_hires.png from https://github.com/JuliaLang/julialang.github.com/tree/master/images

   - Julia Computing
     - jcl-logo.png from http://juliacomputing.com/
     - creators.png from http://venturebeat.com/2015/05/18/why-the-creators-of-the-julia-programming-language-just-launched-a-startup/

   - books.png
     - https://www.packtpub.com/application-development/getting-started-julia-programming/
     - https://pragprog.com/book/7lang/seven-more-languages-in-seven-weeks
     - http://www.manning.com/voncsefalvay/
  """

  md"""
  # (╯°□°）╯  ┻━┻
  """

  md"""
  # end
  """
]

end


text:
using CircGeometry, Test

# colors
julia_blue = :royalblue
julia_green = :forestgreen
julia_red = :brown3
julia_purple = :mediumorchid3
colors = [julia_red,julia_purple,julia_green]

vf = 0.2     #vf = 0.5 for real log (0.5 sometimes times-out with 80 bodies)
n_bodies = 80
material = MaterialParameters(vf,n_bodies)

radius = 0.65
center = Point(0.0,0.0)
outline = OutlineCircle(radius,center)

between_buffer = 15.0
ps1 = generate_porous_structure(outline,material,between_buffer;log=true)
ps2 = generate_porous_structure(outline,material,between_buffer;log=true)
ps3 = generate_porous_structure(outline,material,between_buffer;log=true)

# translate
CircGeometry.translate!(ps1,-0.8,0.0)
CircGeometry.translate!(ps2,0.8,0.0)
CircGeometry.translate!(ps3,0.0,1.2)

# save
save_image("julia_logo.svg",[ps1,ps2,ps3],colors;alpha=1.0)

rm("julia_logo.svg")

max_stars_repo_name: pseastham/CircGeometry

max_issues_repo_name: pseastham/CircGeometry


text:
#File for the Assumption Free K MC2 algorithm (KMeans)
module CustomKMeans

using Distributions
using StatsBase
using Clustering
export KMeansInducingPoints

#Return K inducing points from X, m being the number of Markov iterations for the seeding
function KMeansInducingPoints(X,K,m)
  C = (KmeansSeed(X,K,m))'
  kmeans!(X',C)
  return C'
end
#Fast and efficient seeding for KMeans
function KmeansSeed(X,K,m) #X is the data, K the number of centers wanted, m the number of Markov iterations
  N = size(X,1)
  #Preprocessing, sample first random center
  init = StatsBase.sample(1:N,1)
  C = zeros(K,size(X,2))
  C[1,:] = X[init,:]
  q = zeros(N)
  for i in 1:N
    q[i] = 0.5*norm(X[i,:]-C[1])^2
  end
  sumq = sum(q)
  q = WeightVec(q/sumq + 1.0/(2*N),1)
  uniform = Distributions.Uniform(0,1)
  for i in 2:K
    x = X[StatsBase.sample(1:N,q,1),:] #weighted sampling,
    mindist = mindistance(x,C,i-1)
    for j in 2:m
      y = X[StatsBase.sample(q),:] #weighted sampling
      dist = mindistance(y,C,i-1)
      if (dist/mindist > rand(uniform))
        x = y;  mindist = dist
      end
    end
    C[i,:]=x
  end
  return C
end

#Compue the minimum distance
function mindistance(x,C,K) #Point to look for, collection of centers, number of centers computed
  mindist = Inf
  for i in 1:K
    mindist = min(norm(x-C[i])^2,mindist)
  end
  return mindist
end

end #module


text:
import Logging

@testset "Logger configuration" begin
    logs = let
        logger = TestLogger(min_level=Debug)
        with_logger(logger) do
            @debug "a"
            configure_logging(min_level=Info)
            @debug "a"
            @info  "b"
        end
        logger.logs
    end
    @test length(logs) == 2
    @test occursin((Debug, "a"), logs[1])
    @test occursin((Info , "b"), logs[2])

    # Same test as above, but with global logger
    old_logger = global_logger()
    logs = let
        logger = TestLogger(min_level=Debug)
        global_logger(logger)
        @debug "a"
        configure_logging(min_level=Info)
        @debug "a"
        @info  "b"
        logger.logs
    end
    global_logger(old_logger)

    @test length(logs) == 2
    @test occursin((Debug, "a"), logs[1])
    @test occursin((Info , "b"), logs[2])
end

@testset "configure_logging implementations" begin
    lg1 = ConsoleLogger(stderr, Info)
    lg2 = configure_logging(lg1, min_level=Warn)
    @test lg2.min_level == Warn
    lg3 = configure_logging(lg1, min_level=:warn)
    @test lg3.min_level == Warn

    lg1 = Logging.ConsoleLogger(stderr, Info)
    lg2 = configure_logging(lg1, min_level=Warn)
    @test lg2.min_level == Warn
    lg3 = configure_logging(lg1, min_level=:warn)
    @test lg3.min_level == Warn
end

@testset "disable_logging with parse_level" begin
    # Test utility: Log once at each standard level
    function log_each_level()
        @debug "a"
        @info  "b"
        @warn  "c"
        @error "d"
    end

    disable_logging("Info")
    @test_logs (Warn, "c") (Error, "d")  log_each_level()

    disable_logging(BelowMinLevel)
end


text:
import Base: close
export PeakWriter, close, writePeak, sortPeaks, PeakWriter_unfused, WritePeak_unfused

mutable struct PeakWriter
    #a output stream object to write to
    outstream
    contigs::ReferenceContigs
    cur_ref::Int64
    peakid::Int64
    lastpos::Int64
end

function PeakWriter(output_stream, contigs)
    pw = PeakWriter(output_stream, contigs, 1, 1, 0)
    pw
end

close(pw::PeakWriter) = close(pw.outstream)

function writePeak(pw::PeakWriter, binSize::Int64, binPosStart::Int64, binPosEnd::Int64, pval, fold; prescision=2)

    # Some assertion to prevent misuse.
    @assert pw.lastpos<binPosStart

    # Get the starting position and ending position
    startPos = binSize*(binPosStart-1)+1
    endPos = binSize*binPosEnd

    # If your starting position goes over the size of current chromosome.
    while startPos > pw.contigs.offsets[pw.cur_ref]+pw.contigs.sizes[pw.cur_ref]
        pw.cur_ref += 1
    end

    # Current chrom name
    chr = pw.contigs.names[pw.cur_ref]
    # Get the current position
    startPos = startPos-pw.contigs.offsets[pw.cur_ref]
    endPos = endPos-pw.contigs.offsets[pw.cur_ref]
    # Name peak
    peakname = "peak_$(pw.peakid)"
    # Calculate score (cutoff fold at 1000)
    score = minimum([1000, fold])

    # Write it to file
    output = "$(chr)\t$(startPos)\t$(endPos)\t$(peakname)\t$(round(score,digits=prescision))\t.\t$(round(fold,digits=prescision))\t$(round(pval,digits=prescision))\t-1\t-1"
    println(pw.outstream, output)

    # Update some data
    pw.peakid += 1
    pw.lastpos = binPosEnd

    output
end

#########################################################################
# Preprocesses peak array into list of peaks for peak writer to take in #
#########################################################################
function sortPeaks(pvals, folds, th::Float64)
    flag = false
    pval_list = []
    fold_list = []
    startPos = 0
    peaks = []

    # parse through pvals
    for i in 1:length(pvals)
        if pvals[i] > th

            if flag == false
                # record starting position
                startPos = i
            end

            # set flag
            flag = true
            # record pvals
            push!(pval_list, pvals[i])
            push!(fold_list, folds[i])
        elseif pvals[i] < th && flag
            # set off flag
            flag = false
            # determine max pvals
            maxpval = maximum(pval_list)
            maxfold = maximum(fold_list)
            # reset lists
            pval_list = []
            fold_list = []
            # get ending pos
            endPos = i-1
            # record peaks
            # println([startPos, endPos, maxpval, maxfold])
            push!(peaks, (Int(startPos), Int(endPos), maxpval, maxfold))
        end
    end

    # finish off a peak
    if flag
        maxpval = maximum(pval_list)
        maxfold = maximum(fold_list)
        push!(peaks, (Int(startPos), length(pvals), maxpval, maxfold))
    end

    peaks
end

mutable struct PeakWriter_unfused
    #a output stream object to write to
    Outstream
    contigs::ReferenceContigs
    cur_ref::Int64
    id::Int64
end

function PeakWriter_unfused(output_stream, contigs)
    sw = PeakWriter_unfused(output_stream, contigs, 1, 1)
    sw
end

function WritePeak_unfused(sw::PeakWriter_unfused, binSize::Int64, binPosStart::Int64, binPosEnd::Int64, pval, fold; prescision=3)
    startPos = binSize*(binPosStart-1)+1
    endPos = binSize*binPosEnd

    while startPos > sw.contigs.offsets[sw.cur_ref]+sw.contigs.sizes[sw.cur_ref]
        sw.cur_ref += 1
    end

    chr = sw.contigs.names[sw.cur_ref]
    startPos = startPos-sw.contigs.offsets[sw.cur_ref]
    endPos = endPos-sw.contigs.offsets[sw.cur_ref]
    peakname = "peak_$(sw.id)"
    sw.id += 1
    score = fold

    output = "$(chr)\t$(startPos)\t$(endPos)\t$(peakname)\t$(round(score,digits=prescision))\t.\t$(round(fold,digits=prescision))\t$(round(pval,digits=prescision))\t-1\t-1"

    println(sw.Outstream, output)
end


text:
using Test
using Plots
using LaTeXStrings

@testset "2sls" begin
    include("2sls.jl")
    spacdf[isnan.(spacdf)] .= 0
    @test sum(abs2.(cdf - spacdf)) / length(cdf) < 0.001
end
@testset "es" begin
    include("es.jl")
    @test sum(abs2.(ccdf - ccdfspa)) / length(ccdf) < 0.0001
    @test sum(abs2.(pm - pmspa)) / length(pm) < 0.1
end

@static if ~haskey(ENV, "CI")
    pgfplotsx()
    spacdf[abs.(cdf .- spacdf) .> .1] .= NaN # filter errors due to singularity around the mean
    p1 = plot(xvec, cdf', color=collect(1:length(nuvec))', lw=1.25, legend=:topleft, labels=permutedims("\$\\nu =".*string.(nuvec).*"\$"))
    plot!(
        xvec,
        spacdf',
        color=collect(1:length(nuvec))',
        lw=1.25,
        labels="",
        ls=:dash,
        title="Exact (solid) and approximate (dashes) distributions",
        xlabel=L"$x$",
        ylabel=L"\mathrm{pr}(\hat{\beta}_{2SLS})< x)"
        )
    p2 = plot(ccdf, pm, labels="", lw=1.25)
    plot!(ccdfspa,
        pmspa,
        ls=:dash,
        lw=1.25,
        labels="",
        xlabel=L"VaR level $\alpha$",
        ylabel=L"ES_L^{(\alpha)}",
        title="Exact (solid) and approximate (dashes) expected shortfall"
        )

    plot(p1, p2, layout=@layout[a{1.25w, 1.25h} b{1.25w, 1.25h}])
    savefig("figure1.pdf")
end


text:
module BayesianNetworks

##############################################################################
##
## Dependencies
##
##############################################################################

using Graphs

include("utils.jl")
include("probability_distribution.jl")
include("probability_density_distribution.jl")
include("bayesian_node.jl")
include("bayesian_edge.jl")
include("cpd.jl")
include("bayesian_network.jl")
include("show.jl")

##############################################################################
##
## Exports
##
##############################################################################

export

       # Utils
       normalize_factor,
       normpdf,
       group_combinations,
       joint_probability_distributions,

       # Bayesian Network
       BayesianNetwork,
       add_node!,
       add_edge!,
       nodes,
       num_nodes,
       edges,
       num_edges,
       cpds,
       num_cpds,
       checked_cpds,
       num_checked_cpds,
       set_checked_cpd!,
       find_node_by_symbol,
       symbols_in_network,
       node_in_network,
       nodes_in_network,
       add_probability!,
       legal_configuration,
       check_requirements,
       has_cached_result,

       PDistribution,
       UnknownPDistribution,
       ProbabilityDistribution,
       ProbabilityDensityDistribution,
       states,
       pdfs,
       probabilities,
       CPT,
       BayesianNode,
       DBayesianNode,
       CBayesianNode,
       probability,
       BayesianEdge,
       node_index,
       CPD,
       P,
       valid_pdfs,
       isequal,
       query

end


text:
using Colors

Base.mean{C<:Color}(c1::C,c2::C) = convert(C, (convert(XYZ, c1) + convert(XYZ, c2))/2)  # colorimetrically "proper"

# RGB linear vector space averaging
function meanrgb{C<:AbstractRGB}(c1::C, c2::C)
    C((red(c1)+red(c2))/2,
      (green(c1)+green(c2))/2,
      (blue(c1)+blue(c2))/2)
end

n = 10
colors = distinguishable_colors(n, Array(RGB{Float32},0))
m = 50 # number of repeated pixels in each block
grid = Array(eltype(colors), n*m, n*m)
z = zero(RGB{Float32})
for i = 1:n
    for j = i:n   # i == j is a sanity check
        prop = mean(colors[i], colors[j])
        vs = meanrgb(colors[i], colors[j])
        for l = 1:m, k = 1:m
            grid[(i-1)*m+k,(j-1)*m+l] = grid[(j-1)*m+k,(i-1)*m+l] = k > l ? prop : vs
        end
    end
end
# Put a boundary around the diagonal elements
for i = 1:n
    for k = 1:m
        grid[(i-1)*m+k,(i-1)*m+1] = grid[(i-1)*m+k,(i-1)*m+m] = z
        grid[(i-1)*m+1,(i-1)*m+k] = grid[(i-1)*m+m,(i-1)*m+k] = z
    end
end

using Images

gridu8 = map(Clamp(RGB{U8}), grid)
img = Image(gridu8, spatialorder=["x","y"])
imwrite(img, "comparison.png")


text:
workspace()
using BenchmarkTools
using IntervalArithmetic
using EAGOParametricInterval

# sets options for contractor
opt1 = Any[50    #Int64: Number of iterations
       1.0E-30 #Float64: Tolerance for equality of
       1.0E-10 #Float64: Add Interval(1E-8,1E8) to add to M[i,i] when
              #         processing extended interval division.
      ]
################################################################################
#########################  Test Problem #1 #####################################
################################################################################

P1 = [Interval(5.0,7.0),Interval(5.0,7.0)]
Z1 = [Interval(-1.5, 0.0),Interval(0.0, 0.5)]
h1(z,p) = [z[1]^2+z[2]^2+p[1]*z[1]+4;
           z[1]+p[2]*z[2]]
hj1(z,p) = [(2*z[1]+p[1]) (2*z[2]);
              1              p[2]]

Eflag = false
Iflag = false
eDflag = false
newtonGS1 = PI_NewtonGS(Z1,P1,hj1,h1,opt1,Eflag,Iflag,eDflag)
Eflag = false
Iflag = false
krawczykCW1 = PI_KrawczykCW(Z1,P1,hj1,h1,opt1,Eflag,Iflag)

function h1!(hh,z,p)
       hh[1] = z[1]^2+z[2]^2+p[1]*z[1]+4
       hh[2] = z[1]+p[2]*z[2]
end
function hj1!(hh,z,p)
       hh[1,1] = 2*z[1]+p[1]
       hh[1,2] = 2*z[2]
       hh[2,1] = 1
       hh[2,2] = p[2]
end

Eflag = false
Iflag = false
eDflag = false
InnewtonGS1 = PIn_NewtonGS(Z1,P1,hj1!,h1!,opt1,Eflag,Iflag,eDflag)
Eflag = false
Iflag = false
InkrawczykCW1 = PIn_KrawczykCW(Z1,P1,hj1!,h1!,opt1,Eflag,Iflag)

################################################################################
#########################  Test Problem #2 #####################################
################################################################################

function LBD_func(i)
          YL = [-30, -30, -30, 1800, 900]
          return YL[i]
end
function UBD_func(i)
          YU = [30, 30, 30, 2200, 1100]
          return YU[i]
end

P2 = [Interval(1800,2200),Interval(900,1100)]
Z2 = [Interval(-30, 30),Interval(-30, 30),Interval(-30, 30)]
function h2(x,p)
       return [(3.25-x[1])/p[1]-x[3]; x[1]/p[2]-x[3]; x[2]-(x[1]^2)/(1+x[1]^2)]
end
function hj2(x,p)
       return [-1/p[1] 0 -1;
       1/p[2] 0 -1;
       -2*x[1]/(1+x[1]^2)^2 1 0]
end

Eflag = false
Iflag = false
eDflag = false
newtonGS2 = PI_NewtonGS(Z2,P2,hj2,h2,opt1,Eflag,Iflag,eDflag)
P2 = [Interval(1800,2200),Interval(900,1100)]
Z2 = [Interval(-30, 30),Interval(-30, 30),Interval(-30, 30)]
Eflag = false
Iflag = false
krawczykCW2 = PI_KrawczykCW(Z2,P2,hj2,h2,opt1,Eflag,Iflag)

function h2!(hout,x,p)
       hout[1] = (3.25-x[1])/p[1]-x[3]
       hout[2] = x[1]/p[2]-x[3]
       hout[3] = x[2]-(x[1]^2)/(1+x[1]^2)
end
function hj2!(hout,x,p)
       hout[1,1] = -1/p[1]
       hout[1,2] = 0
       hout[1,3] = -1
       hout[2,1] = 1/p[2]
       hout[2,2] = 0
       hout[2,3] = -1
       hout[3,1] = -2*x[1]/(1+x[1]^2)^2
       hout[3,2] = 1
       hout[3,3] = 0
end
P2 = [Interval(1800,2200),Interval(900,1100)]
Z2 = [Interval(-30, 30),Interval(-30, 30),Interval(-30, 30)]
Eflag = false
Iflag = false
eDflag = false
InnewtonGS2 = PIn_NewtonGS(Z2,P2,hj2!,h2!,opt1,Eflag,Iflag,eDflag)
P2 = [Interval(1800,2200),Interval(900,1100)]
Z2 = [Interval(-30, 30),Interval(-30, 30),Interval(-30, 30)]
Eflag = false
Iflag = false
InkrawczykCW2 = PIn_KrawczykCW(Z2,P2,hj2!,h2!,opt1,Eflag,Iflag)


text:
module Clapeyron
using StaticArrays
using LinearAlgebra
using NLopt, NLSolvers,Roots
#using NLSolve
using  DiffResults, ForwardDiff
include("solvers/Solvers.jl")
using .Solvers
using Combinatorics
using NamedArrays
import LogExpFunctions
using DataStructures: DefaultDict
include("constants.jl")
include("models/basetools.jl")
include("utils/ClapeyronParam.jl")

include("utils/macros.jl")
using CSV, Tables
include("utils/database.jl")
include("utils/misc.jl")

include("models/combiningrules.jl")

include("models/eos.jl")
include("utils/visualisation.jl")

include("models/eos/ideal/BasicIdeal.jl")
include("models/eos/ideal/MonomerIdeal.jl")
include("models/eos/ideal/ReidIdeal.jl")
include("models/eos/ideal/WalkerIdeal.jl")

include("models/eos/SAFT/PCSAFT/PCSAFT.jl")
include("models/eos/SAFT/PCSAFT/variants/sPCSAFT.jl")

include("models/eos/SAFT/ogSAFT/ogSAFT.jl")
include("models/eos/SAFT/SAFTVRSW/SAFTVRSW.jl")
include("models/eos/SAFT/LJSAFT/LJSAFT.jl")
include("models/eos/SAFT/softSAFT/softSAFT.jl")
include("models/eos/SAFT/SAFTVRMie/SAFTVRMie.jl")
include("models/eos/SAFT/SAFTVRMie/variants/SAFTVRQMie.jl")
include("models/eos/SAFT/SAFTgammaMie/SAFTgammaMie.jl")

include("models/eos/SAFT/CKSAFT/CKSAFT.jl")
include("models/eos/SAFT/CKSAFT/variants/sCKSAFT.jl")

include("models/eos/SAFT/BACKSAFT/BACKSAFT.jl")

include("models/eos/cubic/vdW.jl")
include("models/eos/cubic/RK.jl")
include("models/eos/cubic/SRK.jl")
include("models/eos/cubic/PR.jl")
include("models/eos/cubic/CPA.jl")
include("models/eos/cubic/equations.jl")

include("models/eos/EmpiricHelmholtz/IAPWS95.jl")
include("models/eos/EmpiricHelmholtz/PropaneRef.jl")
include("models/eos/EmpiricHelmholtz/GERG2008/GERG2008.jl")

include("models/eos/SPUNG/SPUNG.jl")

# include("models/param_structs.jl")
# include("models/ideal_param_structs.jl")
# include("models/model_structs.jl")
# include("models/eos/combining_rules.jl")
# include("models/import_params.jl")
# include("models/import_ideal_params.jl")

# include("models/eos/Ideal/Basic.jl")
# include("models/eos/Ideal/Monomer.jl")
# include("models/eos/Ideal/Walker.jl")
# include("models/eos/Ideal/Reid.jl")

# include("models/eos/SAFT/PCSAFT.jl")
# include("models/eos/SAFT/CKSAFT.jl")
# include("models/eos/SAFT/sCKSAFT.jl")
# include("models/eos/SAFT/BACKSAFT.jl")
# include("models/eos/SAFT/sPCSAFT.jl")
# include("models/eos/SAFT/SAFTVRMie.jl")
# include("models/eos/SAFT/SAFTVRQMie.jl")
# include("models/eos/SAFT/SAFTVRSW.jl")
# include("models/eos/SAFT/ogSAFT.jl")
# include("models/eos/SAFT/LJSAFT.jl")
# include("models/eos/SAFT/softSAFT.jl")
# include("models/eos/SAFT/SAFTgammaMie.jl")

# include("models/eos/Cubic/vdW.jl")
# include("models/eos/Cubic/RK.jl")
# include("models/eos/Cubic/SRK.jl")
# include("models/eos/Cubic/PR.jl")
# include("models/eos/Cubic/CPA.jl")

# include("models/eos/eos.jl")

# include("models/system.jl")
# include("models/system2.jl")

# export system,System
# export eos,ideal

using Unitful

include("methods/methods.jl")

# export get_volume, get_sat_pure, get_bubble_pressure, get_crit_pure, get_enthalpy_vap, get_pressure, get_entropy, get_chemical_potential, get_internal_energy, get_enthalpy, get_Gibbs_free_energy, get_Helmholtz_free_energy, get_isochoric_heat_capacity, get_isobaric_heat_capacity, get_isothermal_compressibility, get_isentropic_compressibility, get_speed_of_sound, get_isobaric_expansivity, get_Joule_Thomson_coefficient, get_second_virial_coeff
# export eos
# export create_z
end # module


text:
# """
# 	$(SIGNATURES)

# Read a file by [x, gpa]. GPA is in the columns.
# Return the moments for the gpa marginal.
# """
# function read_by_x_gpa(ds :: DataSettings, target :: Symbol)
#     df = read_delim_file_to_df(target);
#     colV = col_headers(:gpa, 1 : n_gpa(ds));
#     m = all_row_to_vector(df, colV);
#     @argcheck length(m) == n_gpa(ds)
#     return m :: Vector{Double}
# end

# Entry rate by GPA
function frac_enter_by_gpa(ds :: DataSettings)
    # Read COL totals b/c file is transposed
    load_fct =
        mt -> read_col_totals(raw_entry_gpa_parental(ds; momentType = mt));
    m, ses, cnts = choice_prob_from_xy(load_fct);
    # m, ses, cnts = choice_prob_from_col_total(ds, :fracEnter_gV);
    # m = read_col_totals(raw_file_path(ds, :fracEnter_gV));
    @assert all(m .> 0.0)  &&  all(m .< 1.0)
    @assert length(m) == n_gpa(ds)

    # cnts = read_col_totals(raw_file_path(ds, :fracEnter_gV; momentType = :count));
    @assert all(cnts .> 100)
    # ses = (m .* (1.0 .- m) ./ cnts) .^ 0.5;
    return m, ses, cnts
end

"""
	$(SIGNATURES)

Mean time to graduation by GPA.
"""
function grad_rate_by_gpa(ds :: DataSettings)
    load_fct =
        mt -> read_col_totals(raw_grad_rate_qual_gpa(ds; momentType = mt));
    m, ses, cnts = choice_prob_from_xy(load_fct);
    # m, ses, cnts = choice_prob_from_col_total(ds, :fracGrad_gV);
    # m = read_col_totals(raw_file_path(ds, :fracGrad_gV));
    @assert all(m .> 0.0)  &&  all(m .< 1.0)
    @assert length(m) == n_gpa(ds)
    @assert all(cnts .> 100)
    return m, ses, cnts
end

# Mean time to drop by quality / gpa
# Contains one 0 cell that needs to be interpolated. 4Y only.
function time_to_drop_4y_by_gpa(ds :: DataSettings)
    load_fct =
        mt -> read_col_totals(ds, :timeToDrop4y_qgM, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert all(m .> 1.5)  &&  all(m .< 4.0)
    @assert length(m) == n_gpa(ds)
    @assert all(cnts .> 30)  "Low counts: $cnts";
    return m, ses, cnts
end

"""
	$(SIGNATURES)

Mean time to graduation by GPA. 4Y colleges only.
"""
function time_to_grad_4y_by_gpa(ds :: DataSettings)
    load_fct =
        mt -> read_col_totals(ds, :timeToGrad4y_qgM, mt);
    m, ses, cnts = load_mean_ses_counts(load_fct);
    @assert all(m .> 3.0)  &&  all(m .< 7.0)
    @assert length(m) == n_gpa(ds)
    @assert all(cnts .>= 30)  "Low counts: $cnts";
    return m, ses, cnts
end

"""
	$(SIGNATURES)

Work hours, year 1, by GPA.
"""
function work_hours_by_gpa(ds :: DataSettings)
    m = read_col_totals(raw_work_hours_qual_gpa(ds, ds.workTimeYear));
    stdV = read_col_totals(
        raw_work_hours_qual_gpa(ds, ds.workTimeYear; momentType = :std));
    cnts = read_col_totals(
        raw_work_hours_qual_gpa(ds, ds.workTimeYear; momentType = :count));
    ses = stdV ./ (max.(cnts, 1.0) .^ 0.5);
    cnts = round.(Int, cnts);
    @assert all(m .> 400.0)  &&  all(m .< 1800.0)
    @assert length(m) == n_gpa(ds)
    @assert all(cnts .> 100)
    return m, ses, cnts
end

## ----------  Dropout rates by GPA / year

# Fraction of initial entrants dropping out at end of each year.
# Standard errors are questionable. The `N`s are given as the total number of students in each college in year 1.
# Dropout rates for 2y starters only sum to about 0.9. But 1/3 occur after year 2.
function frac_drop_gpa_year(ds :: DataSettings)
    # target = :fracDrop_gtM;

    T = ds.Tmax;
    outM = zeros(n_colleges(ds), T);
    sesM = zeros(n_colleges(ds), T);
    cntsM = zeros(Int, n_gpa(ds), T);
    for t = 1 : T
        outM[:, t], sesM[:,t], cntsM[:,t] = frac_drop_gpa(ds, t);
    end
    return outM, sesM, cntsM
end

function frac_drop_gpa(ds :: DataSettings, t :: Integer)
    # Data files contains dropout rates by AFQT as col totals.
    load_fct =
        mt -> read_col_totals(raw_frac_drop_qual_gpa(ds, t; momentType = mt));
    m, ses, cnts = choice_prob_from_xy(load_fct);
    @assert check_float_array(m, 0.0, 1.0)
    @assert length(m) == n_gpa(ds)
    return m, ses, cnts
end

# ------------------


text:
# Showing stuff like power hierarchies
export plot_power_hierarchy, load_cultural_dim, plot_cultural_dim, plot_all_cultural_dim
using CairoMakie, ColorSchemes, ElectronDisplay, Colors

noto_sans = assetpath("fonts", "NotoSans-Regular.ttf")
noto_sans_bold = assetpath("fonts", "NotoSans-Bold.ttf")

function load_cultural_dim(countries::AbstractVector)
    path = joinpath(datadir, "geerthofstede-2015-08-16.csv")
    cultural = CSV.File(path) |> DataFrame
    cultures = select(cultural, :country,
               :pdi => :hierarchy,
               :idv => :individualism,
               :mas => :masculinity,
               :uai => :uncertainty,
               :ltowvs => :long_term,
               :ivr => :indulgence)

    selectedcountries = DataFrame(country=countries)
    sample = innerjoin(cultures, selectedcountries, on = :country)
    select!(sample, [:country, :hierarchy, :individualism, :masculinity, :uncertainty, :long_term, :indulgence])
    sort!(sample, [:hierarchy, :individualism])

    sample
end

function plot_power_hierarchy()
    # countries = ["Norway", "Sweden", "Denmark", "Netherlands", "Germany", "Great Britain", "China", "Iran", "Greece", "India", "Arab countries"]
    # countries = ["Norway", "Sweden", "Denmark", "Netherlands", "Germany", "China", "Iran", "India", "Arab countries"]
    # countries = ["Sweden", "Denmark", "Netherlands", "China", "Iran", "India", "Arab countries"]

    countries = ["Norway", "Denmark", "Germany"]
    df = load_cultural_dim(countries)

    plot_cultural_dim(df, :hierarchy)
end

function plot_cultural_dim(df::DataFrame, column::Symbol)
    xs = 1:length(df.country)
    ys = df[!, column]

    fig = Figure(resolution = (1024, 600), font = noto_sans)
    ax = Axis(fig[1, 1], xticks = (xs, df.country), ylabel=uppercasefirst(string(column)))

    barplot!(ax, xs, ys, color=xs)

    fig
end

"Create legend elements of given colors"
function legend_elements(colors::AbstractArray{<:Color})
    elements = [PolyElement(polycolor = color) for color in colors]
end

function plot_all_cultural_dim()
    countries = ["Norway", "Denmark", "Netherlands", "Germany"]
    df = load_cultural_dim(countries)
    plot_all_cultural_dim(df)
end

function plot_all_cultural_dim(df::DataFrame)
    ncountries = nrow(df)
    # Will interpolate between these 4 colors to create enough colors for ncountries
    colors = get(ColorSchemes.Spectral_4, 1:ncountries, (1, ncountries))

    #culturaldim = ["hierarchy", "individualism", "masculinity", "uncertainty"]
    culturaldim = ["hierarchy", "individualism", "masculinity"]

    xs    = repeat(1:length(culturaldim), inner=ncountries)
    order = repeat(1:ncountries, length(culturaldim))
    #ys    = vcat(df.hierarchy, df.individualism, df.masculinity, df.uncertainty)
    ys    = vcat(df.hierarchy, df.individualism, df.masculinity)

    fig = Figure(resolution = (1120, 630), font = noto_sans)
    ax = Axis(fig[1, 1],
            xticks = (unique(xs), culturaldim))

    barplot!(ax, xs, ys,
        strokewidth=1,
        color = colors[order],
        dodge = order)

    elements = legend_elements(colors)
    title = "Countries"

    Legend(fig[1,2], elements, df.country, title)

    fig
end