text:
reshapeVec2Mat(vec::Vector, rows::Int) = reshape(vec, rows, round(Int,length(vec)/rows))
# function reshapeVec2Mat(vec::Vector, rows::Int)
#   M = reshape(vec, rows, round(Int,length(vec)/rows))
#   return ndims(M) < 2 ? (M')' : M
# end

getData(v::Graphs.ExVertex) = v.attributes["data"]

function getVal(v::Graphs.ExVertex)
  return v.attributes["data"].val
  # return v.attributes["val"]
end

# Convenience function to get values for given variable label
function getVal(fgl::FactorGraph, lbl::ASCIIString)
  getVal(dlapi.getvertex(fgl, lbl))
end
function setVal!(v::Graphs.ExVertex, val::Array{Float64,2})
  v.attributes["data"].val = val
  # v.attributes["val"] = val
  nothing
end
function getBWVal(v::Graphs.ExVertex)
  return v.attributes["data"].bw
end
function setBW!(v::Graphs.ExVertex, bw::Array{Float64,2})
  v.attributes["data"].bw = bw
  # v.attributes["bw"] = bw
  nothing
end
function setVal!(v::Graphs.ExVertex, val::Array{Float64,2}, bw::Array{Float64,2})
  setVal!(v,val)
  setBW!(v,bw)
  nothing
end
function setVal!(v::Graphs.ExVertex, val::Array{Float64,2}, bw::Array{Float64,1})
  setVal!(v,val,(bw')')
  nothing
end
function setValKDE!(v::Graphs.ExVertex, val::Array{Float64,2})
  p = kde!(val)
  setVal!(v,val,getBW(p)[:,1]) # TODO -- this can be little faster
  nothing
end

# function setDefaultNodeDataOld!(v::Graphs.ExVertex, initval::Array{Float64,2},
#                               stdev::Array{Float64,2}, dodims::Int64, N::Int64)
#   pN = Union{}
#   if size(initval,2) < N
#     p = kde!(initval,diag(stdev));
#     pN = resample(p,N)
#   else
#     pN = kde!(initval, "lcv")
#   end
#   # v.attributes["initval"] = initval
#   # v.attributes["initstdev"] = stdev
#   # v.attributes["eliminated"] = false
#   v.attributes["BayesNetVert"] = Union{}
#   dims = size(initval,1) # rows indicate dimensions
#   # v.attributes["dims"] = dims
#   # v.attributes["dimIDs"] = round(Int64,linspace(dodims,dodims+dims-1,dims))
#
#   setDefaultNodeDataReplace!(v,initval,stdev,dodims,N)
#   # setVal!(v, getPoints(pN), getBW(pN)[:,1])
#   nothing
# end

function setDefaultNodeData!(v::Graphs.ExVertex, initval::Array{Float64,2},
                              stdev::Array{Float64,2}, dodims::Int64, N::Int64)
  pN = Union{}
  if size(initval,2) < N
    p = kde!(initval,diag(stdev));
    pN = resample(p,N)
  else
    pN = kde!(initval, "lcv")
  end
  dims = size(initval,1) # rows indicate dimensions
  sp = round(Int64,linspace(dodims,dodims+dims-1,dims))
  data = VariableNodeData(initval, stdev, getPoints(pN),
                          (getBW(pN)[:,1]')', Int64[], sp,
                          dims, false, 0, Int64[])

  v.attributes["data"] = data
  # v.attributes["BayesNetVert"] = Union{} # TODO -- remove
  nothing
end

function addNewVarVertInGraph!(fgl::FactorGraph, vert::Graphs.ExVertex, id::Int64, lbl::AbstractString)
  vert.attributes = Graphs.AttributeDict() #fg.v[fg.id]
  vert.attributes["label"] = lbl #fg.v[fg.id]
  fgl.v[id] = vert # TODO -- this is likely not required, but is used in subgraph methods
  fgl.IDs[lbl] = id # fg.id, to help find it
  nothing
end

# Must rething the abstraction here
function addNode!(fg::FactorGraph, lbl, initval=[0.0]', stdev=[1.0]'; N::Int=100)
  fg.id+=1
  vert = ExVertex(fg.id,lbl)
  addNewVarVertInGraph!(fg, vert, fg.id, lbl)
  # dlapi.setupvertgraph!(fg, vert, fg.id, lbl) #fg.v[fg.id]
  dodims = fg.dimID+1
  # TODO -- vert should not loose information here
  setDefaultNodeData!(vert, initval, stdev, dodims, N) #fg.v[fg.id]

  dlapi.addvertex!(fg, vert) #fg.g ##vertr =

  fg.dimID+=size(initval,1) # rows indicate dimensions, move to last dimension
  push!(fg.nodeIDs,fg.id)
  return vert #fg.v[fg.id]
end

# rethink abstraction, maybe closer to CloudGraph use case a better solution
# function addEdge!(g::FGG,n1,n2)
#   edge = dlapi.makeedge(g, n1, n2)
#   dlapi.addedge!(g, edge)
#   # edge = Graphs.make_edge(g, n1, n2)
#   # Graphs.add_edge!(g, edge)
# end

function getVal(vA::Array{Graphs.ExVertex,1})
  len = length(vA)
  vals = Array{Array{Float64,2},1}()
  cols = Array{Int64,1}()
  push!(cols,0)
  rows = Array{Int64,1}()
  for v in vA
      push!(vals, getVal(v))
      c = size(vals[end],2)
      r = size(vals[end],1)
      push!(cols, floor(Int64,c))
      push!(rows, floor(Int64,r))
  end
  cols = cumsum(cols)
  sc = cols[end]
  rw = floor(Int64,rows[1])
  val = Array{Float64,2}(rw, sc)
  for i in 1:(len-1)
      val[:,(cols[i]+1):cols[i+1]] = vals[i]
  end
  val[:,(cols[len]+1):cols[len+1]] = vals[len] # and the last one
  return val
end

# function getDev(v)
#   return v.attributes["stdev"]
# end

function evalFactor(fg::FactorGraph, evalstr::AbstractString)
  return eval(parse(evalstr))
end

function evalFactor(fg::FactorGraph, fct::Graphs.ExVertex)
  return evalFactor(fg, fct.attributes["evalstr"] )
end

function setDefaultFactorNode!(fact::Graphs.ExVertex, f::Union{Pairwise,Singleton})
  # fact.attributes["fnc"] = f
  # fact.attributes["fncargvID"] = Dict{Int,Int}()
  # fact.attributes["eliminated"] = false
  # fact.attributes["potentialused"] = false

  data = FunctionNodeData{typeof(f)}(Int64[], false, false, Int64[], f)
  fact.attributes["data"] = data

  # for graphviz drawing
  fact.attributes["shape"] = "point"  # fg.f[fg.id]
  fact.attributes["width"] = 0.2 #  fg.f[fg.id]

  nothing
end
# function setFncArgIDs!(fact::Graphs.ExVertex, idx::Int64, index::Int64)
#   push!(fact.attributes["data"].fncargvID, index) #[idx]=index
#   #fact.attributes["fncargvID"][idx] = index
#   nothing
# end

function addNewFncVertInGraph!(fgl::FactorGraph, vert::Graphs.ExVertex, id::Int64, lbl::AbstractString)
  vert.attributes = Graphs.AttributeDict() #fg.v[fg.id]
  vert.attributes["label"] = lbl #fg.v[fg.id]
  fgl.f[id] = vert # TODO -- not sure if this is required
  fgl.fIDs[lbl] = id # fg.id
  nothing
end

function addFactor!(fg::FactorGraph, Xi::Array{Graphs.ExVertex,1},f::Union{Pairwise,Singleton})
  namestring = ""
  for vert in Xi #f.Xi
    namestring = string(namestring,vert.attributes["label"])
  end
  fg.id+=1

  newvert = ExVertex(fg.id,namestring)
  addNewFncVertInGraph!(fg, newvert, fg.id, namestring)
  setDefaultFactorNode!(newvert, f) # fg.f[fg.id]
  push!(fg.factorIDs,fg.id)

  for vert in Xi
    push!(newvert.attributes["data"].fncargvID, vert.index)
  end

  # TODO -- multiple accesses to DB with this method, must refactor!
  newvert = dlapi.addvertex!(fg, newvert)  # used to be two be three lines up ##fg.g
  # idx=0
  for vert in Xi #f.Xi
    dlapi.makeaddedge!(fg, vert, newvert)
    # edge = dlapi.makeedge(fg.g, vert, newvert)
    # dlapi.addedge!(fg.g, edge)
    ## addEdge!(fg.g, vert, newvert) #fg.f[fg.id]
    ## idx+=1
    ## setFncArgIDs!(newvert, idx, vert.index) #fg.f[fg.id]

    #-- push!(newvert.attributes["data"].fncargvID, vert.index)
  end

  return newvert # fg.f[fg.id]
end

function emptyFactorGraph()
    fg = FactorGraph(Graphs.inclist(Graphs.ExVertex,is_directed=false),
                     Graphs.inclist(Graphs.ExVertex,is_directed=true),
                     Dict{Int,Graphs.ExVertex}(),
                     Dict{Int,Graphs.ExVertex}(),
                     Dict{AbstractString,Int}(),
                     Dict{AbstractString,Int}(),
                     0,
                     [],
                     [],
                     Dict{Int,Graphs.ExVertex}(),
                     0,
                     0,
                     nothing,
                     Dict{Int64,Int64}() )
    return fg
end

function prtslperr(s)
  println(s)
  sleep(0.1)
  error(s)
end

# for computing the Bayes Net-----------------------------------------------------
function getEliminationOrder(fg::FactorGraph; ordering::Symbol=:qr)
    s = fg.nodeIDs
    sf = fg.factorIDs
    A=convert(Array{Int},adjacency_matrix(fg.g))[sf,s] # TODO -- order seems brittle
    p = Int[]
    if ordering==:chol
      p = cholfact(A'A,:U,Val{true})[:p] #,pivot=true
    elseif ordering==:qr
      q,r,p = qr(A,Val{true})
    else
      prtslperr("getEliminationOrder -- cannot do the requested ordering $(ordering)")
    end

    # we need the IDs associated with the Graphs.jl and our Type fg
    return fg.nodeIDs[p]
end

# lets create all the vertices first and then deal with the elimination variables thereafter
function addBayesNetVerts!(fg::FactorGraph, elimOrder::Array{Int64,1})
  for p in elimOrder
    vert = dlapi.getvertex(fg, p)
    if vert.attributes["data"].BayesNetVertID == 0   #fg.v[p].
      fg.bnid+=1
      vert.attributes["data"].BayesNetVertID = p # fg.v[p] ##fg.bnverts[p]
      dlapi.updatevertex!(fg, vert)
      # fg.bnverts[p] = Graphs.add_vertex!(fg.bn, ExVertex(fg.bnid,string("BayesNet",fg.bnid)))
      # fg.bnverts[p].attributes["label"] = fg.v[p].attributes["label"]
    else
      println("addBayesNetVerts -- something is very wrong, should not have a Bayes net vertex")
    end
  end
end

function addConditional!(fg::FactorGraph, vertID::Int64, lbl, Si)
  bnv = dlapi.getvertex(fg, vertID) #fg.v[vertID]
  bnvd = bnv.attributes["data"]
  bnvd.separator = Si
  # bnv.attributes["separator"] = Si # TODO -- remove
  # bnvert = bnv.attributes["BayesNetVert"] # TODO -- remove
  for s in Si
    push!(bnvd.BayesNetOutVertIDs, s)
    # addEdge!(fg.bn, fg.v[s].attributes["BayesNetVert"], bnvert) # TODO -- remove
  end
  dlapi.updatevertex!(fg, bnv)
  nothing
end

function addChainRuleMarginal!(fg::FactorGraph, Si)
  lbls = ASCIIString[]
  genmarg = GenericMarginal()
  Xi = Graphs.ExVertex[] #genmarg.Xi = Graphs.ExVertex[]
  for s in Si
    # push!(lbls,fg.v[s].attributes["label"])
    push!(Xi, dlapi.getvertex(fg, s))  # push!(Xi,fg.v[s])
  end
  println("adding marginal to")
  for x in Xi @show x.index end
  # addFactor!(fg, marg, lbls)
  addFactor!(fg, Xi, genmarg)
  nothing
end

# TODO -- Cannot have any CloudGraph calls at this level, must refactor
function rmVarFromMarg(fgl::FactorGraph, fromvert::Graphs.ExVertex, gm::Array{Graphs.ExVertex,1})
  for m in gm
    # get all out edges
    # get neighbors
    for n in dlapi.outneighbors(fgl, m)
      if n.index == fromvert.index
        alleids = m.attributes["data"].edgeIDs
        # println("consider marginal, edge of interest between $(m.index) and $(n.index)")
        i = 0
        for id in alleids
          # println("rmVarFromMarg -- at id=$(id)")
          i+=1
          edge = dlapi.getedge(fgl, id)
          if edge != nothing # hack to avoid dictionary use case
            if edge.SourceVertex.exVertexId == m.index || edge.DestVertex.exVertexId == m.index
              warn("removing edge $(edge.neo4jEdgeId), between $(m.index) and $(n.index)")
              dlapi.deleteedge!(fgl, edge)
              m.attributes["data"].edgeIDs = alleids[[collect(1:(i-1));collect((i+1):length(alleids))]]
              dlapi.updatevertex!(fgl, m)
            end
          end
        end
      end
    end
    # if 0 edges, delete the marginal
    if length(dlapi.outneighbors(fgl, m)) <= 1
      warn("removing vertex id=$(m.index)")
      dlapi.deletevertex!(fgl,m)
    end
  end
  nothing
end

function buildBayesNet!(fg::FactorGraph, p::Array{Int,1})
    addBayesNetVerts!(fg, p)
    for v in p
      println()
      println("Eliminating $(v)")
      println("===============")
      println()
      # which variable are we eliminating
      #fg.v[v].attributes["label"]

      # all factors adjacent to this variable
      fi = Int64[]
      Si = Int64[]
      gm = ExVertex[]
      # TODO -- optimize outneighbor calls like this
      vert = dlapi.getvertex(fg,v)
      for fct in dlapi.outneighbors(fg, vert) #out_neighbors(dlapi.getvertex(fg,v),fg.g) # (fg.v[v]
        if (fct.attributes["data"].eliminated != true) #if (fct.attributes["eliminated"] != true)
          push!(fi, fct.index)
          for sepNode in dlapi.outneighbors(fg, fct) #out_neighbors(fct,fg.g)
            if (sepNode.index != v && length(findin(sepNode.index,Si)) == 0)
              push!(Si,sepNode.index)
            end
          end
          fct.attributes["data"].eliminated = true #fct.attributes["eliminated"] = true
          dlapi.updatevertex!(fg, fct)
        end

        if typeof(fct.attributes["data"].fnc) == GenericMarginal
          push!(gm, fct)
        end
      end

      if v != p[end]
        addConditional!(fg, v, "", Si)
        # not yet inserting the new prior p(Si) back into the factor graph
      end

      tuv = dlapi.getvertex(fg,v)
      tuv.attributes["data"].eliminated = true # fg.v[v].
      dlapi.updatevertex!(fg, tuv)
      # dlapi.getvertex(fg,v).attributes["data"].eliminated = true # fg.v[v].
      ## fg.v[v].attributes["eliminated"] = true

      # TODO -- remove links from current vertex to any marginals
      rmVarFromMarg(fg, vert, gm)

      #add marginal on remaining variables... ? f(xyz) = f(x | yz) f(yz)
      # new function between all Si
      addChainRuleMarginal!(fg, Si)

      # readline(STDIN)
    end
    nothing
end

# some plotting functions on the factor graph
function stackVertXY(fg::FactorGraph, lbl::ASCIIString)
    v = dlapi.getvertex(fg,lbl) # v = fg.v[fg.IDs[lbl]]
    vals = getVal(v)
    X=vec(vals[1,:])
    Y=vec(vals[2,:])
    # X=vec(v.attributes["val"][1,:])
    # Y=vec(v.attributes["val"][2,:])
    return X,Y
end

function getKDE(v::Graphs.ExVertex)
  #return kde!(getVal(v), "lcv") # TODO - getBW(v)
  return kde!(getVal(v), getBWVal(v)[:,1])
end

function getVertKDE(fgl::FactorGraph, id::Int64)
  v = dlapi.getvertex(fgl,id)  # v = fgl.v[fgl.IDs[lbl]]
  return getKDE(v)
end
function getVertKDE(fgl::FactorGraph, lbl::ASCIIString)
  v = dlapi.getvertex(fgl,lbl)  # v = fgl.v[fgl.IDs[lbl]]
  return getKDE(v)
end

function drawCopyFG(fgl::FactorGraph)
  fgd = deepcopy(fgl)
  for v in fgd.v
    delete!(v[2].attributes,"data")
    # delete!(v[2].attributes,"initstdev")
    # delete!(v[2].attributes,"dimIDs")
    # delete!(v[2].attributes,"eliminated")
    # delete!(v[2].attributes,"val")
    # delete!(v[2].attributes,"bw")
    # delete!(v[2].attributes,"BayesNetVert")
    # delete!(v[2].attributes,"initval")
    # delete!(v[2].attributes,"dims")
  end
  for v in fgd.f
    delete!(v[2].attributes,"data")
    # delete!(v[2].attributes,"fncargvID")
    # delete!(v[2].attributes,"eliminated")
    # delete!(v[2].attributes,"fnc")
    # delete!(v[2].attributes,"potentialused")
  end
  return fgd
end

function writeGraphPdf(fgl::FactorGraph)
  fgd = drawCopyFG(fgl)
  println("Writing factor graph file")
  fid = open("fg.dot","w+")
  write(fid,to_dot(fgd.g))
  close(fid)
  run(`dot fg.dot -Tpdf -o fg.pdf`)
  nothing
end

function appendUseFcts!(usefcts, lblid::Int, fct::Graphs.ExVertex, fid::Int)
  for tp in usefcts
    if tp[2].label == fct.label
      # println("Skipping repeat of $(fct.label)")
      return nothing
    end
  end
  tpl = (fct.index, fct) #(lblid, fct, fid)
  push!(usefcts, tpl )
  nothing
end

function expandEdgeListNeigh!(fgl::FactorGraph,
                              vertdict::Dict{Int64,Graphs.ExVertex},
                              edgedict::Dict{Int64,Graphs.Edge{Graphs.ExVertex}})
  #asfd
  for vert in vertdict
    for newedge in out_edges(vert[2],fgl.g)
      if !haskey(edgedict, newedge.index)
        edgedict[newedge.index] = newedge
      end
    end
  end

  nothing
end

# dictionary of unique vertices from edgelist
function expandVertexList!(fgl::FactorGraph,
  edgedict::Dict{Int64,Graphs.Edge{Graphs.ExVertex}},
  vertdict::Dict{Int64,Graphs.ExVertex})

  # go through all source and target nodes
  for edge in edgedict
    if !haskey(vertdict, edge[2].source.index)
      vertdict[edge[2].source.index] = edge[2].source
    end
    if !haskey(vertdict, edge[2].target.index)
      vertdict[edge[2].target.index] = edge[2].target
    end
  end
  nothing
end

function edgelist2edgedict(edgelist::Array{Graphs.Edge{Graphs.ExVertex},1})
  edgedict = Dict{Int64,Graphs.Edge{Graphs.ExVertex}}()
  for edge in edgelist
    edgedict[edge.index] = edge
  end
  return edgedict
end

function addVerticesSubgraph(fgl::FactorGraph,
    fgseg::FactorGraph,
    vertdict::Dict{Int64,Graphs.ExVertex})

    for vert in vertdict
      fgseg.g.vertices[vert[1]] = vert[2]
      if haskey(fgl.v,vert[1])
        fgseg.v[vert[1]] = vert[2]
        fgseg.IDs[vert[2].label] = vert[1]

        # add edges going in opposite direction
        elr = Graphs.out_edges(vert[2], fgl.g)
        len = length(elr)
        keeprm = trues(len)
        j = 0
        for i in 1:len
          if !haskey(vertdict, elr[i].target.index) # a function node in set, so keep ref
            keeprm[i] = false
            j+=1
          end
        end
        if j < len
          elridx = elr[1].source.index
          fgseg.g.inclist[elridx] = elr[keeprm]
        end
      elseif haskey(fgl.f, vert[1])
        fgseg.f[vert[1]] = vert[2] # adding element to subgraph
        fgseg.fIDs[vert[2].label] = vert[1]
        # get edges associated with function nodes and push edges onto incidence list
        el = Graphs.out_edges(vert[2], fgl.g)
        elidx = el[1].source.index
        fgseg.g.inclist[elidx] = el # okay because treating function nodes only
        fgseg.g.nedges += length(el)
      else
        error("Unknown type factor graph vertex type, something is wrong")
      end
    end
    nothing
end

# NOTICE, nodeIDs and factorIDs are not brough over by this method yet
# must sort out for incremental updates
function genSubgraph(fgl::FactorGraph,
    vertdict::Dict{Int64,Graphs.ExVertex})
    # edgedict::Dict{Int64,Graphs.Edge{Graphs.ExVertex}},

  fgseg = FactorGraph() # new handle for just a segment of the graph
  fgseg.g = Graphs.inclist(Graphs.ExVertex,is_directed=false)

  fgseg.v = Dict{Int,Graphs.ExVertex}()
  fgseg.f = Dict{Int,Graphs.ExVertex}()
  fgseg.IDs = Dict{AbstractString,Int}()
  fgseg.fIDs = Dict{AbstractString,Int}()

  fgseg.g.vertices = Array{Graphs.ExVertex,1}(length(fgl.g.vertices))
  fgseg.g.inclist = Array{Array{Graphs.Edge{Graphs.ExVertex},1},1}(length(fgl.g.inclist))

  addVerticesSubgraph(fgl, fgseg, vertdict)

  fgseg.id = fgl.id
  fgseg.bnid = fgl.bnid
  fgseg.dimID = fgl.dimID

  return fgseg
end

function getShortestPathNeighbors(fgl::FactorGraph;
    from::Graphs.ExVertex=nothing,
    to::Graphs.ExVertex=nothing,
    neighbors::Int=0 )

  edgelist = shortest_path(fgl.g, ones(num_edges(fgl.g)), from, to)
  vertdict = Dict{Int64,Graphs.ExVertex}()
  edgedict = edgelist2edgedict(edgelist)
  expandVertexList!(fgl, edgedict, vertdict) # grow verts
  for i in 1:neighbors
    # println("Expanding subgraph edge dict")
    expandEdgeListNeigh!(fgl, vertdict, edgedict) # grow edges
    expandVertexList!(fgl, edgedict, vertdict) # grow verts
  end
  return vertdict
end

function subgraphShortestPath(fgl::FactorGraph;
    from::Graphs.ExVertex=nothing,
    to::Graphs.ExVertex=nothing,
    neighbors::Int=0  )

  vertdict = getShortestPathNeighbors(fgl, from=from, to=to, neighbors=neighbors)
  return genSubgraph(fgl, vertdict)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
function subgraphFromVerts(fgl::FactorGraph,
    verts::Dict{Int64,Graphs.ExVertex};
    neighbors::Int=0  )

    allverts = Dict{Int64,Graphs.ExVertex}()
    allkeys = collect(keys(verts))
    len = length(allkeys)
    # union all shortest path combinations in a vertdict
    for i in 1:len, j in (i+1):len
      from = verts[allkeys[i]]
      to = verts[allkeys[j]]
      vertdict = getShortestPathNeighbors(fgl, from=from, to=to, neighbors=neighbors)
      for vert in vertdict
        if !haskey(allverts, vert[1])
          allverts[vert[1]] = vert[2]
        end
      end
    end

  return genSubgraph(fgl, allverts)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
# Using unique index into graph data structure
function subgraphFromVerts(fgl::FactorGraph,
    verts::Array{Int64,1};
    neighbors::Int=0  )

  vertdict = Dict{Int,Graphs.ExVertex}()
  for vert in verts
    vertdict[vert] = fgl.v[vert]
  end

  return subgraphFromVerts(fgl,vertdict,neighbors=neighbors)
end

# explore all shortest paths combinations in verts, add neighbors and reference subgraph
# Using unique index into graph data structure
function subgraphFromVerts(fgl::FactorGraph,
    verts::Array{ASCIIString,1};
    neighbors::Int=0  )

  vertdict = Dict{Int,Graphs.ExVertex}()
  for vert in verts
    id = -1
    if haskey(fgl.IDs, vert)
      id = fgl.IDs[vert]
    else
      error("FactorGraph01 only supports variable node subgraph search at this point")
    end
    vertdict[id] = fgl.v[id]
  end

  return subgraphFromVerts(fgl,vertdict,neighbors=neighbors)
end


text:
function ls_params(ϕ::AbstractVector{T}, c::Bool, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where T
    P = isempty(p) ? 0 : maximum(p)
    Np = length(p)
    Q = isempty(q) ? 0 : maximum(q)
    Nq = length(q)
    ϕ0 = zeros(T, P)
    θ0 = zeros(T, Q)
    μ0 = c ? ϕ[1] : zero(T)
    for (i, l) in enumerate(p)
        ϕ0[l] = ϕ[i + c]
    end
    for (i, l) in enumerate(q)
        θ0[l] = ϕ[i + Np + c]
    end
    β0 = ϕ[Nq + Np + 1 + c:end]
    return μ0, ϕ0, θ0, β0
end

function ls_matrix!(Z::AbstractMatrix{T}, z::AbstractVector{T}, a::AbstractVector{T}, x::AbstractMatrix{T}, c::Bool, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where T
    all(q .!= 0) || throw(ArgumentError("order of ma can't be zero"))
    all(p .!= 0) || throw(ArgumentError("order of ar can't be zero"))
    N = length(z)
    P = isempty(p) ? 0 : maximum(p)
    Np = length(p)
    Q = isempty(q) ? 0 : maximum(q)
    Nq = length(q)
    X = size(x, 2)
    if c
        Z[:,1] .= 1
    end
    for (i, l) in enumerate(p)
        @views Z[:,i + c] = z[P + Q - l + 1:N - l]
    end
    for (i, l) in enumerate(q)
        @views Z[:,i + Np + c] = -a[P + Q - l + 1:N - l]
    end
    if X > 0 @views Z[:, Np + Nq + 1 + c:end] = x[P:N - 1, :] end
    return Z, z[P + Q + 1:N]
end

function ls_matrix(z::AbstractVector{T}, a::AbstractVector{T}, x::AbstractMatrix{T}, c::Bool, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where T
    N = length(z)
    P = isempty(p) ? 0 : maximum(p)
    Np = length(p)
    Q = isempty(q) ? 0 : maximum(q)
    Nq = length(q)
    X = size(x, 2)
    Z = Matrix{T}(undef, N - P - Q, Np + Nq + X + c)
    ls_matrix!(Z, z, a, x, c, p, q)
end

function least_squares(z::AbstractVector{T}, a::AbstractVector{T}, x::AbstractMatrix{T}, c::Bool, p::AbstractVector{<:Integer}, q::AbstractVector{<:Integer}) where T
    Z, z = ls_matrix(z, a, x, c, p, q)
    ϕ = Z \ z
    ε = z - Z * ϕ
    σ2 = dot(ε, ε) / length(ε)
    μ, ϕ, θ, β = ls_params(ϕ, c, p, q)
    return μ, ϕ, θ, β, σ2
end

least_squares(z::AbstractVector{T}, a::AbstractVector{T}, x::AbstractMatrix{T}, c::Bool, p::Integer, q::Integer) where T = least_squares(z, a, x, c, collect(0:p), collect(1:q))

ext: jl


text:
using Documenter, NewGuilia

makedocs(;
    modules=[NewGuilia],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/DarioSarra/NewGuilia.jl/blob/{commit}{path}#L{line}",
    sitename="NewGuilia.jl",
    authors="Dario Sarra",
    assets=String[],
)

deploydocs(;
    repo="github.com/DarioSarra/NewGuilia.jl",
)


text:
include("setup_1d.jl")

println("idwt by filtering (N=",N,"), ", L, " levels")
f(x0, L, wf) = for i = 1:tn; idwt(x0, wf, L); end
f(x0, L, wf);
@time f(x0, L, wf);


text:
using Bokeh

plotfile("hold.html")

x = linspace(-3pi, 3pi)
y1 = cos(x)
plot(x,y1)

hold(true)

y2 = sin(x)./x
plot(x, y2, "r")
showplot()


text:
"""
	Peak(f, dBGain, BWoct, fs=96e3)

Return Biquad of sampling frequency fs as DSP.SecondOrderSections that resembles Lake Peak EQ.

The algorithm is a modified version of the method first presented in:

Orfanidis, JAES 1997 Issue 6: "Digital Parametric Equalizer Design with Prescribed Nyquist-Frequency Gain"
http://www.aes.org/e-lib/browse.cfm?elib=7854

In contrast to the method presented in the paper, this function does not
define linear Δf based on Bandwidth but forces the level to be half of
dBGain at the lower Bandwidth-defining point, e. g. 1 octave down from f when
the Bandwidth is 2 octaves. The resulting Biquads correlated very well with
measured frequency responses from the Lake Controller.
"""
function Peak(f, dBGain, BWoct, fs=96e3)
	if dBGain==0
		return SecondOrderSections(ZeroPoleGain(Float64[], Float64[], 1.))
	end
	G_0 = 1 #gain at f=0
	G = 10^(dBGain/20) #gain at peak
	G_b = 10^(dBGain/40) #gain that defines bandwidth
	f_b = f * 2^(-BWoct/2) #lower frequency that defines bandwidth
	#calculate digital frequencies
	omega_0 = f/fs*2pi
	omega_b = f_b/fs*2pi
	#prewarping
	Omega_0 = tan(omega_0/2)
	Omega_b = tan(omega_b/2)
	#Nyquist-Gain from analog model
	#F is frequency in rad
	F_s = fs*pi #Nyquist frequency
	F_0 = f*2pi #peak frequency
	F_1 = F_0*2^(-BWoct/2) #lower frequency that defines bandwidth
	F_2 = F_0*2^(+BWoct/2) #upper frequency that defines bandwidth
	DeltaF = F_2-F_1
	A_analog = sqrt((G_b^2-G_0^2)/(G^2-G_b^2))*DeltaF
	B_analog = G * A_analog
	G_1 = sqrt((G_0*(F_s^2-F_0^2)^2+B_analog^2*F_s^2)/((F_s^2-F_0^2)^2+A_analog^2*F_s^2))
	#analog parameters to be bilinearly transformed, imported from wxMaxima
	A = sqrt(-(((G_0^2-G^2)*G_b^2+(G^2-G_0^2)*G_1^2)*Omega_b^4+((2*G^2-2*G_0^2)*sqrt((G_1^2-G^2)/(G_0^2-G^2))*G_b^2+(2*G^2*G_0^2-2*G^4)*sqrt((G_1^2-G^2)/(G_0^2-G^2))+(2*G_0^2-2*G^2)*G_1^2-2*G^2*G_0^2+2*G^4)*Omega_0^2*Omega_b^2+((G_1^2-G^2)*G_b^2-G_0^2*G_1^2+G^2*G_0^2)*Omega_0^4)/((G_0^2-G^2)*G_b^2-G^2*G_0^2+G^4))/Omega_b
	B = sqrt(-(((G^2*G_0^2-G^4)*G_b^2+(G^4-G^2*G_0^2)*G_1^2)*Omega_b^4+(((2*G^2*G_0-2*G_0^3)*G_1*sqrt((G_1^2-G^2)/(G_0^2-G^2))+(2*G_0^2-2*G^2)*G_1^2-2*G^2*G_0^2+2*G^4)*G_b^2+(2*G^2*G_0^3-2*G^4*G_0)*G_1*sqrt((G_1^2-G^2)/(G_0^2-G^2)))*Omega_0^2*Omega_b^2+((G^2*G_1^2-G^4)*G_b^2-G^2*G_0^2*G_1^2+G^4*G_0^2)*Omega_0^4)/((G_0^2-G^2)*G_b^2-G^2*G_0^2+G^4))/Omega_b
	W2 = ((G_1^2-G^2)/(G_0^2-G^2))^(1/2)*Omega_0^2 #W2 = W^2
	#bilinear transform
	b = [
		(G_1+G_0*W2+B),
		-2*(G_1-G_0*W2),
		(G_1-B+G_0*W2)]./(1+W2+A)
	a = [
		1,
		-2*(1-W2)/(1+W2+A),
		(1+W2-A)/(1+W2+A)]
	PeakFilter = SecondOrderSections(Biquad(b[1], b[2], b[3], a[2], a[3]))
	return PeakFilter
end


text:
using LogarithmicNumbers, Test

_approx(x,y) = (x ≈ y) || (isnan(x) && isnan(y))

# use these to check for type stability
_exp(args...) = @inferred exp(args...)
_log(args...) = @inferred log(args...)
_add(args...) = @inferred +(args...)
_sub(args...) = @inferred -(args...)
_mul(args...) = @inferred *(args...)
_div(args...) = @inferred /(args...)
_pow(args...) = @inferred ^(args...)
_float(args...) = @inferred float(args...)
_inv(args...) = @inferred inv(args...)
_prod(args...) = @inferred prod(args...)
_sum(args...) = @inferred sum(args...)

# sample values
vals = (-Inf, -100, -100.0, -2, -2.0, -1, -1.0, -0.5, 0, 0.0, 0.5, 1, 1.0, 2, 2.0, 100, 100.0, Inf)

# sample vectors
vecs = (
	Float64[x for x in vals if x isa Float64],
	Float64[x for x in vals if x isa Float64 && !isinf(x)],
	Float64[x for x in vals if x isa Float64 && !iszero(x)],
	Float64[x for x in vals if x isa Float64 && !isinf(x) && x ≥ 0],
	Int[x for x in vals if x isa Int],
	Int[x for x in vals if x isa Int && x ≥ 0],
)

atypes = (ULogarithmic, Logarithmic, CLogarithmic)

@testset "LogarithmicNumbers" begin

	@testset "types" begin
	    @test @isdefined ULogarithmic
	    @test @isdefined Logarithmic
	    @test @isdefined CLogarithmic
	end

	@testset "exp" begin
		for A in atypes, x in vals
			@test _exp(A, x) isa A{typeof(x)}
		end
	end

	@testset "convert" begin
	    for A in atypes, x in vals
	    	if x < 0 && A == ULogarithmic
	    		@test_throws DomainError A(x)
	    	else
	    		@test (@inferred A(x)) isa A
	    	end
	    end
	end

    @testset "log" begin
    	for A in atypes, x in vals
    		y = A==CLogarithmic ? Complex(x) : x
			@test _log(_exp(A, x)) === y
    	end
    end

	@testset "float" begin
	    for A in atypes, x in vals
	    	_float(_exp(A, x)) ≈ exp(x) || @info("float1",A,x)
	    	@test _float(_exp(A, x)) ≈ exp(x)
	    	if x < 0 && A == ULogarithmic
	    		nothing
	    	elseif x==-Inf && A == CLogarithmic
	    		@test AbstractFloat(A(x)) == x
	    	else
	    		@test _float(A(x)) ≈ x
	    	end
	    end
	end

    @testset "mul" begin
    	for A in atypes, x in vals, y in vals
    		@test _approx(_float(_mul(_exp(A, x), _exp(A, y))), exp(x) * exp(y))
    	end
    end

    @testset "div" begin
    	for A in atypes, x in vals, y in vals
    		@test _approx(_float(_div(_exp(A, x), _exp(A, y))), exp(x) / exp(y))
    	end
    end

    @testset "pow" begin
		for A in atypes, x in vals, n in (-2,-1,0,1,2,-1.1,0.0,2.3)
			if isinf(x) && iszero(n)
				# should be 1.0 but is NaN
				@test_broken _approx(_float(_pow(_exp(A, x), n)), exp(x)^n)
				@test isnan(_float(_pow(_exp(A, x), n)))
			else
			    @test _approx(_float(_pow(_exp(A, x), n)), exp(x)^n)
			end
        end
    end

    @testset "inv" begin
        for A in atypes, x in vals
            @test _approx(_float(_inv(_exp(A, x))), inv(exp(x)))
        end
    end

    @testset "add" begin
        for A in atypes, x in vals, y in vals
    		@test _approx(_float(_add(_exp(A, x), _exp(A, y))), exp(x)+exp(y))
        end
    end

    @testset "sub" begin
        for A in atypes, x in vals, y in vals
        	if A==ULogarithmic && x<y
	    		@test_throws DomainError _float(_sub(_exp(A, x), _exp(A, y)))
	    	elseif A==CLogarithmic && y==Inf && x<y
	    		# problem with complex logarithmics
	    		@test_broken _approx(_float(_sub(_exp(A, x), _exp(A, y))), exp(x)-exp(y))
	    	else
	    		@test _approx(_float(_sub(_exp(A, x), _exp(A, y))), exp(x)-exp(y))
	    	end
        end
    end

    @testset "prod" begin
        for A in atypes, xs in vecs
        	_approx(_float(_prod(map(x->_exp(A,x), xs))), prod(exp.(xs))) || @info("prod",A,xs)
        	@test _approx(_float(_prod(map(x->_exp(A,x), xs))), prod(exp.(xs)))
        end
    end

    @testset "sum" begin
        for A in atypes, xs in vecs
        	if eltype(xs) <: AbstractFloat
	            @test _approx(_float(_sum(map(x->exp(A,x),xs))), sum(exp.(xs)))
	        else
	        	# sum is not type-stable because typeof(xs[1]+xs[2]) != typeof(xs[1]).
	        	# hence this test is the same as above but without the stability check.
	        	# we don't promise type stability unless the base type is a float, so
	        	# this isn't a broken test.
	        	@test _approx(_float(sum(map(x->exp(A,x),xs))), sum(exp.(xs)))
	        end
        end
    end

end


text:
module SquarenessBenchmark

import Squareness
using Squareness: issquare_with_isqrt
using BenchmarkTools

issquare_with_padic_newton(n::Signed) = n < 0 ? false : Squareness.issquare_with_padic_newton(n % Unsigned)
# issquare_with_padic_newton(n::Unsigned) = Squareness.issquare_with_padic_newton(n)
issquare_with_padic_newton(n::BigInt) = n < 0 ? false : Squareness.issquare_with_padic_newton(n)

function _issquare_with_simple_newton(x::I) where {I <: Integer}
    x == 0 && return true
    xtz = trailing_zeros(x)
    isodd(xtz) && return false
    x >>= unsigned(xtz)
    o_xrt::I = x
    xrt::I = (x + 0x1) >> 0x1
    xrt == 0 && return false  # overflow
    while xrt < o_xrt
        o_xrt = xrt
        xrt = (xrt + x ÷ xrt) >> 0x1
    end
    xrt * xrt == x
end
issquare_with_simple_newton(x::Signed) = x < 0 ? false : _issquare_with_simple_newton(x % Unsigned)
# issquare_with_simple_newton(x::Unsigned) = _issquare_with_simple_newton(x)
issquare_with_simple_newton(x::BigInt) = x < 0 ? false : _issquare_with_simple_newton(x)

println("## Int64 ≥ 0\n")

rand_int64() = rand(Int64) & typemax(Int64)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton(n) setup=(n=rand_int64()))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt(n) setup=(n=rand_int64()))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton(n) setup=(n=rand_int64()))
println("\n```\n\n")

println("## Int64 ≥ 0 × 100\n")

rand_int64(N) = rand(Int64, N) .& typemax(Int64)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton.(n) setup=(n=rand_int64(100)))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt.(n) setup=(n=rand_int64(100)))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton.(n) setup=(n=rand_int64(100)))
println("\n```\n\n")

println("## Int32 ≥ 0\n")

rand_int32() = rand(Int32) & typemax(Int32)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton(n) setup=(n=rand_int32()))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt(n) setup=(n=rand_int32()))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton(n) setup=(n=rand_int32()))
println("\n```\n\n")

println("## Int32 ≥ 0 × 100\n")

rand_int32(N) = rand(Int32, N) .& typemax(Int32)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton.(n) setup=(n=rand_int32(100)))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt.(n) setup=(n=rand_int32(100)))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton.(n) setup=(n=rand_int32(100)))
println("\n```\n\n")

println("## Int16 ≥ 0\n")

rand_int16() = rand(Int16) & typemax(Int16)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton(n) setup=(n=rand_int16()))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt(n) setup=(n=rand_int16()))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton(n) setup=(n=rand_int16()))
println("\n```\n\n")

println("## Int16 ≥ 0 × 100\n")

rand_int16(N) = rand(Int16, N) .& typemax(Int16)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton.(n) setup=(n=rand_int16(100)))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt.(n) setup=(n=rand_int16(100)))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton.(n) setup=(n=rand_int16(100)))
println("\n```\n\n")

println("## Int8 ≥ 0\n")

rand_int8() = rand(Int8) & typemax(Int8)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton(n) setup=(n=rand_int8()))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt(n) setup=(n=rand_int8()))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton(n) setup=(n=rand_int8()))
println("\n```\n\n")

println("## Int8 ≥ 0 × 100\n")

rand_int8(N) = rand(Int8, N) .& typemax(Int8)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton.(n) setup=(n=rand_int8(100)))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt.(n) setup=(n=rand_int8(100)))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton.(n) setup=(n=rand_int8(100)))
println("\n```\n\n")

println("## Int128 ≥ 0\n")

rand_int128() = rand(Int128) & typemax(Int128)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton(n) setup=(n=rand_int128()))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt(n) setup=(n=rand_int128()))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton(n) setup=(n=rand_int128()))
println("\n```\n\n")

println("## Int128 ≥ 0 × 100\n")

rand_int128(N) = rand(Int128, N) .& typemax(Int128)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton.(n) setup=(n=rand_int128(100)))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt.(n) setup=(n=rand_int128(100)))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton.(n) setup=(n=rand_int128(100)))
println("\n```\n\n")

println("## BigInt ≥ 0\n")

rand_mpz() = (big(rand(UInt128)) << rand(0:128)) + rand(UInt128)

println("### P-adic Newton Method\n")
println("```")
display(@benchmark issquare_with_padic_newton(n) setup=(n=rand_mpz()))
println("\n```\n\n")

println("### With `isqrt()`\n")
println("```")
display(@benchmark issquare_with_isqrt(n) setup=(n=rand_mpz()))
println("\n```\n\n")

println("### Simple Newton Method\n")
println("```")
display(@benchmark issquare_with_simple_newton(n) setup=(n=rand_mpz()))
println("\n```\n\n")

end


text:
#GracePlot constant literals
#-------------------------------------------------------------------------------

#==Julia Symbol => Grace constants
===============================================================================#

#A constant litteral in grace...
type GraceConstLitteral
	#Basically just a string, but will not be surrounded with quotes when sent...
	s::AbstractString
end

const graceconstmap = Dict{Symbol, GraceConstLitteral}(
	#Booleans:
	:on  => GraceConstLitteral("ON"),
	:off => GraceConstLitteral("OFF"),
	:TRUE  => GraceConstLitteral("TRUE"), #Julia does not like :true/:false
	:FALSE => GraceConstLitteral("FALSE"),

	#Axis scales:
	:lin        => GraceConstLitteral("NORMAL"),
	:log        => GraceConstLitteral("LOGARITHMIC"),
	:reciprocal => GraceConstLitteral("RECIPROCAL"),

	#Common:
	:none       => GraceConstLitteral("0"), #Linestyle, ..

	#Line styles:
	:solid       => GraceConstLitteral("1"),
	:dot         => GraceConstLitteral("2"),
	:dash        => GraceConstLitteral("3"),
	:ldash       => GraceConstLitteral("4"),
	:dotdash     => GraceConstLitteral("5"),
	:dotldash    => GraceConstLitteral("6"),
	:dotdotdash  => GraceConstLitteral("7"),
	:dotdashdash => GraceConstLitteral("8"),

	#Symobls (glyphs):
	:circle    => GraceConstLitteral("1"),
	:o         => GraceConstLitteral("1"),
	:square    => GraceConstLitteral("2"),
	:diamond   => GraceConstLitteral("3"),
	:uarrow    => GraceConstLitteral("4"),
	:larrow    => GraceConstLitteral("5"),
	:darrow    => GraceConstLitteral("6"),
	:rarrow    => GraceConstLitteral("7"),
	:cross     => GraceConstLitteral("8"),
	:+         => GraceConstLitteral("8"),
	:diagcross => GraceConstLitteral("9"),
	:x         => GraceConstLitteral("9"),
	:star      => GraceConstLitteral("10"),
	:*         => GraceConstLitteral("10"),
	:char      => GraceConstLitteral("11"),

	#Location type (relative to what coordinates):
	:world     => GraceConstLitteral("world"),
	:view      => GraceConstLitteral("view"),

	#Text justification
	:centercenter => GraceConstLitteral("14"),
)

#Last line


text:
module Test_MEcondest

using LinearAlgebra
using LinearMaps
using MatrixEquations
using Test

function check_ctranspose(op::LinearMaps.LinearMap{T}) where {T <: Union{AbstractFloat, Complex}}
   (m, n) = size(op)
   x = rand(T,n)
   y = rand(T,m)
   yAx = dot(y, op * x)
   xAty = dot(x, op' * y)
   ε = eps(real(eltype(op)))
   return abs(yAx - conj(xAty)) < (abs(yAx) + ε) * ε^(1 / 3)
end

@testset "Testing Lyapunov and Sylvester operators" begin

n = 10; m = 5;
reltol = sqrt(eps(10.));
sc = 0.1;

@testset "Continuous Lyapunov operators" begin

#n = 3
ar = rand(n,n);
cr = rand(n,m);
cc = cr+im*rand(n,m);
cr = cr*cr';
as,  = schur(ar);
ac = ar+im*rand(n,n);
cc = cc*cc';
acs,  = schur(ac);

Tcr = lyapop(ar);
Tcrinv = invlyapop(ar);
Tcrs = lyapop(as);
Tcrsinv = invlyapop(as);
Tcrsym = lyapop(ar,her=true);
Tcrsyminv = invlyapop(ar,her=true);
Tcrssym = lyapop(as,her=true);
Tcrssyminv = invlyapop(as,her=true);

Tcc = lyapop(ac);
Tccinv = invlyapop(ac);
Tccs = lyapop(acs);
Tccsinv = invlyapop(acs);
Tccsym = lyapop(ac,her=true);
Tccsyminv = invlyapop(ac,her=true);
Tccssym = lyapop(acs,her=true);
Tccssyminv = invlyapop(acs,her=true);
Π = trmatop(n);
@test issymmetric(Π) && ishermitian(Π)
@test size(Tcr) == size(Tcr') == size(Tcrinv)

@test check_ctranspose(Tcr) &&
      check_ctranspose(Tcrinv) &&
      check_ctranspose(Tcrs) &&
      check_ctranspose(Tcrsinv) &&
      #check_ctranspose(Tcrsym) &&
      #check_ctranspose(Tcrsyminv) &&
      #check_ctranspose(Tcrssyminv) &&
      #check_ctranspose(Tcrssyminv) &&
      check_ctranspose(Tcc) &&
      check_ctranspose(Tccinv) &&
      check_ctranspose(Tccs) &&
      check_ctranspose(Tccsinv)
      #check_ctranspose(Tccsym) &&
      #check_ctranspose(Tccsyminv) &&
      #check_ctranspose(Tccssyminv) &&
      #check_ctranspose(Tccssyminv)

@test opnorm1(invlyapop([ 0 0;0 0])) == Inf

try
   T1 = invlyapop([1 0;0 1])
   T1 = invlyapop(convert(Matrix{Float32},as))
   T = invlyapop(as);
   T*rand(n*n);
   transpose(T)*rand(n*n);
   adjoint(T)*rand(n*n);
   T*rand(Float32,n*n);
   transpose(T)*rand(Float32,n*n);
   adjoint(T)*rand(Float32,n*n);
   T*ones(Int,n*n);
   transpose(T)*ones(Int,n*n);
   adjoint(T)*ones(Int,n*n);
   T*ones(Rational{Int},n*n);
   T1 = invlyapop(convert(Matrix{Complex{Float32}},acs));
   T = invlyapop(acs);
   T*rand(n*n);
   T*complex(rand(n*n));
   transpose(T)*rand(n*n);
   transpose(T)*complex(rand(n*n));
   adjoint(T)*rand(n*n);
   adjoint(T)*complex(rand(n*n));
   lyapop(schur(ar))
   lyapop(2)
   @test true
catch
   @test false
end

@time x = Tcrinv*cr[:];
@test norm(Tcr*x-cr[:])/norm(x[:]) < reltol

@time x = Tcrinv*cc[:];
@test norm(Tcr*x-cc[:])/norm(x[:]) < reltol

@time x = Tcrsyminv*triu2vec(cr);
@test norm(Tcrsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tcrinv)*cr[:];
@test norm(transpose(Tcr)*x-cr[:])/norm(x[:]) < reltol

@time x = transpose(Tcrinv)*cc[:];
@test norm(transpose(Tcr)*x-cc[:])/norm(x[:]) < reltol

@time x = Tcrinv'*cc[:];
@test norm(Tcr'*x-cc[:])/norm(x[:]) < reltol

@time x = transpose(Tcrsyminv)*triu2vec(cr);
@test norm(transpose(Tcrsym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tcrsinv*cr[:];
@test norm(Tcrs*x-cr[:])/norm(x[:]) < reltol

@time x = Tcrsinv*cc[:];
@test norm(Tcrs*x-cc[:])/norm(x[:]) < reltol

@time x = transpose(Tcrsinv)*cr[:];
@test norm(transpose(Tcrs)*x-cr[:])/norm(x[:]) < reltol

@time x = transpose(Tcrsinv)*cc[:];
@test norm(transpose(Tcrs)*x-cc[:])/norm(x[:]) < reltol

@time x = Tcrsinv'*cr[:];
@test norm(Tcrs'*x-cr[:])/norm(x[:]) < reltol

@time x = Tcrsinv'*cc[:];
@test norm(Tcrs'*x-cc[:])/norm(x[:]) < reltol

@time x = Tcrssyminv*triu2vec(cr);
@test norm(Tcrssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tcrssyminv)*triu2vec(cr);
@test norm(transpose(Tcrssym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tccinv*cc[:];
@test norm(Tcc*x-cc[:])/norm(x[:]) < reltol

@time x = Tccsyminv*triu2vec(cr);
@test norm(Tccsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = Tccinv'*cc[:];
@test norm(Tcc'*y-cc[:])/norm(y[:]) < reltol

@time x = Tccsyminv'*triu2vec(cr);
@test norm(Tccsym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tccsinv*cc[:];
@test norm(Tccs*x-cc[:])/norm(x[:]) < reltol

@time x = Tccssyminv*triu2vec(cr);
@test norm(Tccssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = transpose(Tcrinv)*cr[:];
@test norm(transpose(Tcr)*y-cr[:])/norm(y[:]) < reltol

@time y = transpose(Tcrsinv)*cr[:];
@test norm(transpose(Tcrs)*y-cr[:])/norm(y[:]) < reltol

@test abs(norm(Matrix(Tcr)'-Matrix(transpose(Tcr)))) < reltol
@test abs(norm(Matrix(Tcc)'-Matrix(Tcc'))) < reltol
@test norm(Matrix(Tcrinv)*Matrix(Tcr)-I) < reltol
@test norm(Matrix(Tcrsinv)*Matrix(Tcrs)-I) < reltol
@test norm(Matrix(Tccinv)*Matrix(Tcc)-I) < reltol
@test norm(Matrix(Tccsinv)*Matrix(Tccs)-I) < reltol
@test opnorm1est(Π*Tcr-Tcr*Π) < reltol
@test opnorm1est(Π*Tcrinv-Tcrinv*Π) < reltol
@test opnorm1est(Π*Tcrsinv-Tcrsinv*Π) < reltol
@test opnorm(Matrix(Tcr),1) ≈ opnorm1(Tcr)
@test opnorm1(Tcrsinv-invlyapop(as')') == 0
@test opnorm1(Tcrsinv-invlyapop(schur(ar))) == 0

@test sc*opnorm1(Tcr) < opnorm1est(Tcr)  &&
      sc*opnorm1(Tcrinv) < opnorm1est(Tcrinv)  &&
      sc*opnorm1(Tcrs) < opnorm1est(Tcrs)  &&
      sc*opnorm1(Tcrsinv) < opnorm1est(Tcrsinv)  &&
      sc*opnorm1(Tcc) < opnorm1est(Tcc)  &&
      sc*opnorm1(Tccinv) < opnorm1est(Tccinv)  &&
      sc*opnorm1(Tccs) < opnorm1est(Tccs)  &&
      sc*opnorm1(Tccsinv) < opnorm1est(Tccsinv)

@test sc*oprcondest(Tcr,Tcrinv) < 1/opnorm1(Tcr)/opnorm1(Tcrinv)  &&
      sc*oprcondest(Tcrs,Tcrsinv) < 1/opnorm1(Tcrs)/opnorm1(Tcrsinv)  &&
      sc*oprcondest(Tcc,Tccinv) < 1/opnorm1(Tcc)/opnorm1(Tccinv)  &&
      sc*oprcondest(Tccs,Tccsinv) < 1/opnorm1(Tccs)/opnorm1(Tccsinv) &&
      oprcondest(Tcr) == oprcondest(Tcr,Tcrinv) == oprcondest(Tcrinv)

@test opsepest(Tcrsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Tcrs)))  &&
      minimum(svdvals(Matrix(Tcrs))) <= sqrt(2)*n*opsepest(Tcrsinv)  &&
      opsepest(invlyapop([0. 1.; 0. 1.])) == 0. &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.]))) == 0. &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.]))) == 0.  &&
      opsepest(invlyapop([0. 1.; 0. 1.])) == 0.  &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.]))) == 0.  &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.]))) == 0.

end

@testset "Continuous generalized Lyapunov operators" begin

ar = rand(n,n)
er = rand(n,n)
cr = rand(n,m)
cc = cr+im*rand(n,m)
cr = cr*cr'
as,es  = schur(ar,er)
ac = ar+im*rand(n,n)
ec = er+im*rand(n,n)
cc = cc*cc'
acs, ecs  = schur(ac,ec)

Tcr = lyapop(ar,er)
Tcrinv = invlyapop(ar,er)
Tcrs = lyapop(as,es)
Tcrsinv = invlyapop(as,es)
Tcrsym = lyapop(ar,er,her=true);
Tcrsyminv = invlyapop(ar,er,her=true);
Tcrssym = lyapop(as,es,her=true);
Tcrssyminv = invlyapop(as,es,her=true);

Tcc = lyapop(ac,ec)
Tccinv = invlyapop(ac,ec)
Tccs = lyapop(acs,ecs)
Tccsinv = invlyapop(acs,ecs)
Tccsym = lyapop(ac,ec,her=true);
Tccsyminv = invlyapop(ac,ec,her=true);
Tccssym = lyapop(acs,ecs,her=true);
Tccssyminv = invlyapop(acs,ecs,her=true);
Π = trmatop(size(ar))
@test size(Tcr) == size(Tcr') == size(Tcrinv)

@test opnorm1est(Π*Tcrinv-Tcrinv*Π) < reltol*norm(ar)*norm(er)

@test check_ctranspose(Tcr) &&
      check_ctranspose(Tcrinv) &&
      check_ctranspose(Tcrs) &&
      check_ctranspose(Tcrsinv) &&
      #check_ctranspose(Tcrsym) &&
      #check_ctranspose(Tcrsyminv) &&
      #check_ctranspose(Tcrssyminv) &&
      #check_ctranspose(Tcrssyminv) &&
      check_ctranspose(Tcc) &&
      check_ctranspose(Tccinv) &&
      check_ctranspose(Tccs) &&
      check_ctranspose(Tccsinv)
      #check_ctranspose(Tccsym) &&
      #check_ctranspose(Tccsyminv) &&
      #check_ctranspose(Tccssyminv) &&
      #check_ctranspose(Tccssyminv)

try
   T1 = invlyapop(convert(Matrix{Complex{Float32}},acs),ecs);
   T = invlyapop(acs,ecs);
   T*rand(n*n);
   transpose(T)*rand(n*n);
   adjoint(T)*rand(n*n);
   T*rand(Float32,n*n);
   transpose(T)*rand(Float32,n*n);
   adjoint(T)*rand(Float32,n*n);
   T*ones(Int,n*n);
   transpose(T)*ones(Int,n*n);
   adjoint(T)*ones(Int,n*n);
   T1 = invlyapop(convert(Matrix{Complex{Float32}},acs),convert(Matrix{Complex{Float32}},ecs));
   T = invlyapop(acs,ecs);
   T*rand(n*n);
   T*complex(rand(n*n));
   transpose(T)*rand(n*n);
   transpose(T)*complex(rand(n*n));
   T'*rand(n*n);
   T'*complex(rand(n*n));
   lyapop(schur(ar,er))
   lyapop(2,3)
   lyapop(rand(Int,2,2),rand(2,2))
   @test true
catch
   @test false
end

@time x = Tcrinv*cr[:];
@test norm(Tcr*x-cr[:])/norm(x[:]) < reltol

@time x = Tcrinv*cc[:];
@test norm(Tcr*x-cc[:])/norm(x[:]) < reltol

@time x = Tcrsyminv*triu2vec(cr);
@test norm(Tcrsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tcrsyminv*triu2vec(cc);
@test norm(Tcrsym*x-triu2vec(cc))/norm(x[:]) < reltol

@time x = transpose(Tcrinv)*cr[:];
@test norm(transpose(Tcr)*x-cr[:])/norm(x[:]) < reltol

@time x = transpose(Tcrinv)*cc[:];
@test norm(transpose(Tcr)*x-cc[:])/norm(x[:]) < reltol

@time x = adjoint(Tcrinv)*cc[:];
@test norm(adjoint(Tcr)*x-cc[:])/norm(x[:]) < reltol

@time x = transpose(Tcrsyminv)*triu2vec(cr);
@test norm(transpose(Tcrsym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tcrsyminv)*triu2vec(cc);
@test norm(transpose(Tcrsym)*x-triu2vec(cc))/norm(x[:]) < reltol

@time x = Tcrsinv*cr[:];
@test norm(Tcrs*x-cr[:])/norm(x[:]) < reltol

@time x = transpose(Tcrsinv)*cr[:];
@test norm(transpose(Tcrs)*x-cr[:])/norm(x[:]) < reltol

@time x = Tcrssyminv*triu2vec(cr);
@test norm(Tcrssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tcrssyminv)*triu2vec(cr);
@test norm(transpose(Tcrssym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = adjoint(Tcrssyminv)*triu2vec(cr);
@test norm(adjoint(Tcrssym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tccinv*cc[:];
@test norm(Tcc*x-cc[:])/norm(x[:]) < reltol

@time x = Tccsyminv*triu2vec(cr);
@test norm(Tccsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = Tccinv'*cc[:];
@test norm(Tcc'*y-cc[:])/norm(y[:]) < reltol

@time x = Tccsyminv'*triu2vec(cr);
@test norm(Tccsym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tccsinv*cc[:];
@test norm(Tccs*x-cc[:])/norm(x[:]) < reltol

@time x = Tccssyminv*triu2vec(cr);
@test norm(Tccssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = Tcrinv'*cr[:];
@test norm(Tcr'*y-cr[:])/norm(y[:]) < reltol

@time y = Tcrsinv'*cr[:];
@test norm(Tcrs'*y-cr[:])/norm(y[:]) < reltol

@test abs(norm(Matrix(Tcr)'-Matrix(transpose(Tcr)))) < reltol
@test abs(norm(Matrix(Tcc)'-Matrix(Tcc'))) < reltol
@test norm(Matrix(Tcrinv)*Matrix(Tcr)-I) < reltol
@test norm(Matrix(Tcrsinv)*Matrix(Tcrs)-I) < reltol
@test norm(Matrix(Tccinv)*Matrix(Tcc)-I) < reltol
@test norm(Matrix(Tccsinv)*Matrix(Tccs)-I) < reltol
@test opnorm1est(Π*Tcr-Tcr*Π) < reltol
@test opnorm1est(Π*Tcrinv-Tcrinv*Π) < reltol*norm(ar)*norm(er)
@test opnorm1est(Π*Tcrsinv-Tcrsinv*Π) < reltol*norm(ar)*norm(er)
@test opnorm(Matrix(Tcr),1) ≈ opnorm1(Tcr)
@test opnorm1(Tcrsinv-invlyapop(as',es')') == 0
@test opnorm1(Tcrsinv-invlyapop(schur(ar,er))) == 0

@test sc*opnorm1(Tcr) < opnorm1est(Tcr)  &&
      sc*opnorm1(Tcrinv) < opnorm1est(Tcrinv)  &&
      sc*opnorm1(Tcrs) < opnorm1est(Tcrs)  &&
      sc*opnorm1(Tcrsinv) < opnorm1est(Tcrsinv)  &&
      sc*opnorm1(Tcc) < opnorm1est(Tcc)  &&
      sc*opnorm1(Tccinv) < opnorm1est(Tccinv)  &&
      sc*opnorm1(Tccs) < opnorm1est(Tccs)  &&
      sc*opnorm1(Tccsinv) < opnorm1est(Tccsinv)

@test sc*oprcondest(Tcr,Tcrinv) < 1/opnorm1(Tcr)/opnorm1(Tcrinv)  &&
      sc*oprcondest(Tcrs,Tcrsinv) < 1/opnorm1(Tcrs)/opnorm1(Tcrsinv)  &&
      sc*oprcondest(Tcc,Tccinv) < 1/opnorm1(Tcc)/opnorm1(Tccinv)  &&
      sc*oprcondest(Tccs,Tccsinv) < 1/opnorm1(Tccs)/opnorm1(Tccsinv) &&
      oprcondest(Tcr) == oprcondest(Tcr,Tcrinv) == oprcondest(Tcrinv)

@test opsepest(Tcrsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Tcrs)))  &&
      minimum(svdvals(Matrix(Tcrs))) <= sqrt(2)*n*opsepest(Tcrsinv)  &&
      opsepest(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.])) == 0. &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.]))) == 0. &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.]))) == 0.  &&
      opsepest(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.])) == 0.  &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.]))) == 0.  &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.]))) == 0.

end

@testset "Discrete Lyapunov operators" begin

ar = rand(n,n);
cr = rand(n,m);
cc = cr+im*rand(n,m);
cr = cr*cr';
as,  = schur(ar);
ac = ar+im*rand(n,n);
cc = cc*cc';
acs,  = schur(ac);

Tdr = lyapop(ar,disc=true);
Tdrinv = invlyapop(ar,disc=true);
Tdrs = lyapop(as,disc=true);
Tdrsinv = invlyapop(as,disc=true);
Tdrsym = lyapop(ar,disc=true,her=true);
Tdrsyminv = invlyapop(ar,disc=true,her=true);
Tdrssym = lyapop(as,disc=true,her=true);
Tdrssyminv = invlyapop(as,disc=true,her=true);

Tdc = lyapop(ac,disc=true);
Tdcinv = invlyapop(ac,disc=true);
Tdcs = lyapop(acs,disc=true);
Tdcsinv = invlyapop(acs,disc=true);
Tdcsym = lyapop(ac,disc=true,her=true);
Tdcsyminv = invlyapop(ac,disc=true,her=true);
Tdcssym = lyapop(acs,disc=true,her=true);
Tdcssyminv = invlyapop(acs,disc=true,her=true);
Π = trmatop(ar)

@test check_ctranspose(Tdr) &&
      check_ctranspose(Tdrinv) &&
      check_ctranspose(Tdrs) &&
      check_ctranspose(Tdrsinv) &&
      check_ctranspose(Tdc) &&
      check_ctranspose(Tdcinv) &&
      check_ctranspose(Tdcs) &&
      check_ctranspose(Tdcsinv)

try
   T1 = invlyapop(convert(Matrix{Float32},as))
   T = invlyapop(as,disc=true);
   T*rand(n*n);
   transpose(T)*rand(n*n);
   adjoint(T)*rand(n*n);
   T*rand(Float32,n*n);
   transpose(T)*rand(Float32,n*n);
   adjoint(T)*rand(Float32,n*n);
   T*ones(Int,n*n);
   transpose(T)*ones(Int,n*n);
   adjoint(T)*ones(Int,n*n);
   T1 = invlyapop(convert(Matrix{Complex{Float32}},acs));
   T = invlyapop(acs,disc=true);
    T*rand(n*n);
    T*complex(rand(n*n));
    transpose(T)*rand(n*n);
    transpose(T)*complex(rand(n*n));
    adjoint(T)*rand(n*n);
    adjoint(T)*complex(rand(n*n));
    @test true
catch
    @test false
end

@time x = Tdrinv*cr[:];
@test norm(Tdr*x-cr[:])/norm(x[:]) < reltol

@time x = Tdrinv*cc[:];
@test norm(Tdr*x-cc[:])/norm(x[:]) < reltol

@time x = Tdrsyminv*triu2vec(cr);
@test norm(Tdrsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tdrinv)*cr[:];
@test norm(transpose(Tdr)*x-cr[:])/norm(x[:]) < reltol

@time x = Tdrsyminv'*triu2vec(cr);
@test norm(Tdrsym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdrsyminv'*triu2vec(cc);
@test norm(Tdrsym'*x-triu2vec(cc))/norm(x[:]) < reltol

@time x = Tdrsinv*cr[:]
@test norm(Tdrs*x-cr[:])/norm(x[:]) < reltol

@time x = transpose(Tdrsinv)*cr[:]
@test norm(transpose(Tdrs)*x-cr[:])/norm(x[:]) < reltol

@time x = Tdrssyminv*triu2vec(cr);
@test norm(Tdrssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tdrsyminv)*triu2vec(cr);
@test norm(transpose(Tdrsym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdrssyminv'*triu2vec(cr);
@test norm(Tdrssym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdcinv*cc[:]
@test norm(Tdc*x-cc[:])/norm(x[:]) < reltol

@time x = Tdcsyminv*triu2vec(cr);
@test norm(Tdcsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = Tdcinv'*cc[:]
@test norm(Tdc'*y-cc[:])/norm(y[:]) < reltol

@time x = Tdcsyminv'*triu2vec(cr);
@test norm(Tdcsym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdcsinv*cc[:]
@test norm(Tdcs*x-cc[:])/norm(x[:]) < reltol

@time x = Tdcssyminv*triu2vec(cr);
@test norm(Tdcssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = transpose(Tdrinv)*cr[:]
@test norm(transpose(Tdr)*y-cr[:])/norm(y[:]) < reltol

@time y = transpose(Tdrsinv)*cr[:]
@test norm(transpose(Tdrs)*y-cr[:])/norm(y[:]) < reltol

@test abs(norm(Matrix(Tdr)'-Matrix(transpose(Tdr)))) < reltol
@test abs(norm(Matrix(Tdc)'-Matrix(Tdc'))) < reltol
@test norm(Matrix(Tdrinv)*Matrix(Tdr)-I) < reltol
@test norm(Matrix(Tdrsinv)*Matrix(Tdrs)-I) < reltol
@test norm(Matrix(Tdcinv)*Matrix(Tdc)-I) < reltol
@test norm(Matrix(Tdcsinv)*Matrix(Tdcs)-I) < reltol
@test opnorm1est(Π*Tdr-Tdr*Π) < reltol
@test opnorm1est(Π*Tdrinv-Tdrinv*Π) < reltol
@test opnorm1est(Π*Tdrsinv-Tdrsinv*Π) < reltol

@test sc*opnorm1(Tdr) < opnorm1est(Tdr)  &&
      sc*opnorm1(Tdrinv) < opnorm1est(Tdrinv)  &&
      sc*opnorm1(Tdrs) < opnorm1est(Tdrs)  &&
      sc*opnorm1(Tdrsinv) < opnorm1est(Tdrsinv)  &&
      sc*opnorm1(Tdc) < opnorm1est(Tdc)  &&
      sc*opnorm1(Tdcinv) < opnorm1est(Tdcinv)  &&
      sc*opnorm1(Tdcs) < opnorm1est(Tdcs)  &&
      sc*opnorm1(Tdcsinv) < opnorm1est(Tdcsinv)

@test sc*oprcondest(Tdr,Tdrinv) < 1/opnorm1(Tdr)/opnorm1(Tdrinv)  &&
      sc*oprcondest(Tdrs,Tdrsinv) < 1/opnorm1(Tdrs)/opnorm1(Tdrsinv)  &&
      sc*oprcondest(Tdc,Tdcinv) < 1/opnorm1(Tdc)/opnorm1(Tdcinv)  &&
      sc*oprcondest(Tdcs,Tdcsinv) < 1/opnorm1(Tdcs)/opnorm1(Tdcsinv) &&
      oprcondest(Tdr) == oprcondest(Tdr,Tdrinv) == oprcondest(Tdrinv)

@test opsepest(Tdrsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Tdrs)))  &&
      minimum(svdvals(Matrix(Tdrs))) <= sqrt(2)*n*opsepest(Tdrsinv)  &&
      opsepest(invlyapop([0. 1.; 0. 1.],disc=true)) == 0. &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.],disc=true))) == 0. &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.],disc=true))) == 0.  &&
      opsepest(invlyapop([0. 1.; 0. 1.],disc=true)) == 0.  &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.],disc=true))) == 0.  &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.],disc=true))) == 0.

end

@testset "Discrete generalized Lyapunov operators" begin

ar = rand(n,n)
er = rand(n,n)
cr = rand(n,m)
cc = cr+im*rand(n,m)
cr = cr*cr'
as,es  = schur(ar,er)
ac = ar+im*rand(n,n)
ec = er+im*rand(n,n)
cc = cc*cc'
acs, ecs  = schur(ac,ec)

Tdr = lyapop(ar,er,disc=true)
Tdrinv = invlyapop(ar,er,disc=true)
Tdrs = lyapop(as,es,disc=true)
Tdrsinv = invlyapop(as,es,disc=true)
Tdrsym = lyapop(ar,er,disc=true,her=true);
Tdrsyminv = invlyapop(ar,er,disc=true,her=true);
Tdrssym = lyapop(as,es,disc=true,her=true);
Tdrssyminv = invlyapop(as,es,disc=true,her=true);

Tdc = lyapop(ac,ec,disc=true)
Tdcinv = invlyapop(ac,ec,disc=true)
Tdcs = lyapop(acs,ecs,disc=true)
Tdcsinv = invlyapop(acs,ecs,disc=true)
Tdcsym = lyapop(ac,ec,disc=true,her=true);
Tdcsyminv = invlyapop(ac,ec,disc=true,her=true);
Tdcssym = lyapop(acs,ecs,disc=true,her=true);
Tdcssyminv = invlyapop(acs,ecs,disc=true,her=true);
Π = trmatop(n)

@test check_ctranspose(Tdr) &&
      check_ctranspose(Tdrinv) &&
      check_ctranspose(Tdrs) &&
      check_ctranspose(Tdrsinv) &&
      check_ctranspose(Tdc) &&
      check_ctranspose(Tdcinv) &&
      check_ctranspose(Tdcs) &&
      check_ctranspose(Tdcsinv)

try
   #T1 = invlyapsop(convert(Matrix{Complex{Float32}},acs),ecs);
   T = invlyapop(acs,ecs,disc = true);
   T*rand(n*n);
   transpose(T)*rand(n*n);
   adjoint(T)*rand(n*n);
   T*rand(Float32,n*n);
   transpose(T)*rand(Float32,n*n);
   adjoint(T)*rand(Float32,n*n);
   T*ones(Int,n*n);
   transpose(T)*ones(Int,n*n);
   adjoint(T)*ones(Int,n*n);
   #T1 = invlyapsop(convert(Matrix{Complex{Float32}},acs),convert(Matrix{Complex{Float32}},ecs));
   T = invlyapop(acs,ecs,disc = true);
   T*rand(n*n);
   T*complex(rand(n*n));
   transpose(T)*rand(n*n);
   transpose(T)*complex(rand(n*n));
   T'*rand(n*n);
   T'*complex(rand(n*n));
   @test true
catch
   @test false
end

@time x = Tdrinv*cr[:];
@test norm(Tdr*x-cr[:])/norm(x[:]) < reltol

@time x = Tdrsyminv*triu2vec(cr);
@test norm(Tdrsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = transpose(Tdrinv)*cr[:];
@test norm(transpose(Tdr)*x-cr[:])/norm(x[:]) < reltol

@time x = Tdrsyminv'*triu2vec(cr);
@test norm(Tdrsym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdrsinv*cr[:];
@test norm(Tdrs*x-cr[:])/norm(x[:]) < reltol

@time x = transpose(Tdrsinv)*cr[:]
@test norm(transpose(Tdrs)*x-cr[:])/norm(x[:]) < reltol

@time x = Tdrssyminv*triu2vec(cr);
@test norm(Tdrssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdrssyminv*triu2vec(cc);
@test norm(Tdrssym*x-triu2vec(cc))/norm(x[:]) < reltol

@time x = transpose(Tdrsyminv)*triu2vec(cr);
@test norm(transpose(Tdrsym)*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdrssyminv'*triu2vec(cr);
@test norm(Tdrssym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdrssyminv'*triu2vec(cc);
@test norm(Tdrssym'*x-triu2vec(cc))/norm(x[:]) < reltol

@time x = Tdcinv*cc[:];
@test norm(Tdc*x-cc[:])/norm(x[:]) < reltol

@time x = Tdcsyminv*triu2vec(cr);
@test norm(Tdcsym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = Tdcinv'*cc[:];
@test norm(Tdc'*y-cc[:])/norm(y[:]) < reltol

@time x = Tdcsyminv'*triu2vec(cr);
@test norm(Tdcsym'*x-triu2vec(cr))/norm(x[:]) < reltol

@time x = Tdcsinv*cc[:]
@test norm(Tdcs*x-cc[:])/norm(x[:]) < reltol

@time x = Tdcssyminv*triu2vec(cr);
@test norm(Tdcssym*x-triu2vec(cr))/norm(x[:]) < reltol

@time y = transpose(Tdrinv)*cr[:];
@test norm(transpose(Tdr)*y-cr[:])/norm(y[:]) < reltol

@time y = transpose(Tdrsinv)*cr[:];
@test norm(transpose(Tdrs)*y-cr[:])/norm(y[:]) < reltol

@test abs(norm(Matrix(Tdr)'-Matrix(transpose(Tdr)))) < reltol
@test abs(norm(Matrix(Tdc)'-Matrix(Tdc'))) < reltol
@test norm(Matrix(Tdrinv)*Matrix(Tdr)-I) < reltol
@test norm(Matrix(Tdrsyminv)*Matrix(Tdrsym)-I) < reltol
@test norm(Matrix(Tdrsinv)*Matrix(Tdrs)-I) < reltol
@test norm(Matrix(Tdrssyminv)*Matrix(Tdrssym)-I) < reltol
@test norm(Matrix(Tdcinv)*Matrix(Tdc)-I) < reltol
@test norm(Matrix(Tdcsinv)*Matrix(Tdcs)-I) < reltol
@test opnorm1est(Π*Tdr-Tdr*Π) < reltol
@test opnorm1est(Π*Tdrinv-Tdrinv*Π) < reltol
@test opnorm1est(Π*Tdrsinv-Tdrsinv*Π) < reltol

@test sc*opnorm1(Tdr) < opnorm1est(Tdr)  &&
      sc*opnorm1(Tdrinv) < opnorm1est(Tdrinv)  &&
      sc*opnorm1(Tdrs) < opnorm1est(Tdrs)  &&
      sc*opnorm1(Tdrsinv) < opnorm1est(Tdrsinv)  &&
      sc*opnorm1(Tdc) < opnorm1est(Tdc)  &&
      sc*opnorm1(Tdcinv) < opnorm1est(Tdcinv)  &&
      sc*opnorm1(Tdcs) < opnorm1est(Tdcs)  &&
      sc*opnorm1(Tdcsinv) < opnorm1est(Tdcsinv)

@test sc*oprcondest(Tdr,Tdrinv) < 1/opnorm1(Tdr)/opnorm1(Tdrinv)  &&
      sc*oprcondest(Tdrs,Tdrsinv) < 1/opnorm1(Tdrs)/opnorm1(Tdrsinv)  &&
      sc*oprcondest(Tdc,Tdcinv) < 1/opnorm1(Tdc)/opnorm1(Tdcinv)  &&
      sc*oprcondest(Tdcs,Tdcsinv) < 1/opnorm1(Tdcs)/opnorm1(Tdcsinv) &&
      oprcondest(Tdr) == oprcondest(Tdr,Tdrinv) == oprcondest(Tdrinv)

@test opsepest(Tdrsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Tdrs)))  &&
      minimum(svdvals(Matrix(Tdrs))) <= sqrt(2)*n*opsepest(Tdrsinv)  &&
      opsepest(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true)) == 0. &&
      opsepest(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true)) == 0. &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true))) == 0. &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true))) == 0.  &&
      opsepest(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true)) == 0.  &&
      opsepest(transpose(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true))) == 0.  &&
      opsepest(adjoint(invlyapop([0. 1.; 0. 1.],[1. 1.;0. 1.],disc=true))) == 0.

end

#  continuous and discrete Sylvester equations
@testset "Continuous Sylvester operators" begin

# n = 3; m = 2;
ar = rand(n,n)
br = rand(m,m)
cr = rand(n,m)
as,  = schur(ar)
bs,  = schur(br)
ac = ar+im*rand(n,n)
bc = br+im*rand(m,m)
cc = cr+im*rand(n,m)
acs,  = schur(ac)
bcs,  = schur(bc)

Tcr = sylvop(ar, br)
Tcrinv = invsylvop(ar,br)
Tcrs = sylvop(as, bs)
Tcrsinv = invsylvop(as,bs)

Tcc = sylvop(ac, bc)
Tccinv = invsylvop(ac,bc)
Tccs = sylvop(acs, bcs)
Tccsinv = invsylvop(acs,bcs)

@test size(Tcr) == size(Tcr')

@test check_ctranspose(Tcr) &&
      check_ctranspose(Tcrinv) &&
      check_ctranspose(Tcrs) &&
      check_ctranspose(Tcrsinv) &&
      check_ctranspose(Tcc) &&
      check_ctranspose(Tccinv) &&
      check_ctranspose(Tccs) &&
      check_ctranspose(Tccsinv)

try
    T = invsylvop(as,bs);
    T*rand(n*m);
    transpose(T)*rand(n*m);
    adjoint(T)*rand(n*m);
    T*rand(Float32,n*m);
    transpose(T)*rand(Float32,n*m);
    adjoint(T)*rand(Float32,n*m);
    T*ones(Int,n*m);
    transpose(T)*ones(Int,n*m);
    adjoint(T)*ones(Int,n*m);
    T = invsylvop(acs,bcs);
    T*rand(n*m);
    T*complex(rand(n*m));
    transpose(T)*rand(n*m);
    transpose(T)*complex(rand(n*m));
    T'*rand(n*m);
    T'*complex(rand(n*m));
    sylvop(schur(rand(3,3)),schur(rand(3,3)))
    sylvop(1,im)
    invsylvop(1,im)
    sylvop(rand(Int,2,2),rand(2,2))
    invsylvop(rand(Int,2,2),rand(2,2))
    @test true
catch
    @test false
end

@time x = Tcrinv*cr[:];
@test norm(Tcr*x[:]-cr[:])/norm(x[:]) < reltol

@time x = Tcrsinv*cr[:];
@test norm(Tcrs*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as,bs)*cr[:];
@test norm(sylvop(as, bs)*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as,bs')*cr[:];
@test norm(sylvop(as, bs')*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as',bs)*cr[:];
@test norm(sylvop(as', bs)*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as',bs')*cr[:];
@test norm(sylvop(as', bs')*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as,bs))*cr[:];
@test norm(transpose(sylvop(as, bs))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as,bs'))*cr[:];
@test norm(transpose(sylvop(as, bs'))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as',bs))*cr[:];
@test norm(transpose(sylvop(as', bs))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as',bs'))*cr[:];
@test norm(transpose(sylvop(as', bs'))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as,bs))*cr[:];
@test norm(adjoint(sylvop(as, bs))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as,bs'))*cr[:];
@test norm(adjoint(sylvop(as, bs'))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as',bs))*cr[:];
@test norm(adjoint(sylvop(as', bs))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as',bs'))*cr[:];
@test norm(adjoint(sylvop(as', bs'))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = Tccinv*cc[:];
@test norm(Tcc*x[:]-cc[:])/norm(x[:]) < reltol

@time x = Tccsinv*cc[:];
@test norm(Tccs*x[:]-cc[:])/norm(x[:]) < reltol

@time y = transpose(Tcrinv)*cr[:];
@test norm(transpose(Tcr)*y[:]-cr[:])/norm(y[:]) < reltol

@time y = Tccinv'*cc[:];
@test norm(Tcc'*y[:]-cc[:])/norm(y[:]) < reltol

x = rand(n*m);
@test abs(norm(Matrix(Tcr)'-Matrix(transpose(Tcr)))) < reltol
@test abs(norm(Matrix(Tcc)'-Matrix(Tcc'))) < reltol
@test norm(Matrix(Tcrinv)*Matrix(Tcr)-I) < reltol
@test norm(Matrix(Tcrsinv)*Matrix(Tcrs)-I) < reltol
@test norm(Matrix(Tccinv)*Matrix(Tcc)-I) < reltol
@test norm(Matrix(Tccsinv)*Matrix(Tccs)-I) < reltol
@test norm(Tcrsinv'*x-invsylvop(as',bs')*x) == 0
@test opnorm1(Tcrsinv-invsylvop(schur(ar),schur(br))) == 0

@test sc*opnorm1(Tcr) < opnorm1est(Tcr)  &&
      sc*opnorm1(Tcrinv) < opnorm1est(Tcrinv)  &&
      sc*opnorm1(Tcrs) < opnorm1est(Tcrs)  &&
      sc*opnorm1(Tcrsinv) < opnorm1est(Tcrsinv)  &&
      sc*opnorm1(Tcc) < opnorm1est(Tcc)  &&
      sc*opnorm1(Tccinv) < opnorm1est(Tccinv)  &&
      sc*opnorm1(Tccs) < opnorm1est(Tccs)  &&
      sc*opnorm1(Tccsinv) < opnorm1est(Tccsinv)

@test sc*oprcondest(Tcr,Tcrinv) < 1/opnorm1(Tcr)/opnorm1(Tcrinv)  &&
      sc*oprcondest(Tcrs,Tcrsinv) < 1/opnorm1(Tcrs)/opnorm1(Tcrsinv)  &&
      sc*oprcondest(Tcc,Tccinv) < 1/opnorm1(Tcc)/opnorm1(Tccinv)  &&
      sc*oprcondest(Tccs,Tccsinv) < 1/opnorm1(Tccs)/opnorm1(Tccsinv) &&
      oprcondest(Tcr) == oprcondest(Tcr,Tcrinv) == oprcondest(Tcrinv)

@test opsepest(Tcrsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Tcrs)))  &&
      minimum(svdvals(Matrix(Tcrs))) <= sqrt(2)*n*opsepest(Tcrsinv)  &&
      opsepest(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.])) == 0. &&
      opsepest(transpose(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.]))) == 0. &&
      opsepest(adjoint(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.]))) == 0. &&
      opsepest(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.])) == 0. &&
      opsepest(transpose(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.]))) == 0. &&
      opsepest(adjoint(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.]))) == 0.

end

# discrete Sylvester equations
@testset "Discrete Sylvester operators" begin

ar = rand(n,n)
br = rand(m,m)
cr = rand(n,m)
as,  = schur(ar)
bs,  = schur(br)
ac = ar+im*rand(n,n)
bc = br+im*rand(m,m)
cc = cr+im*rand(n,m)
acs,  = schur(ac)
bcs,  = schur(bc)

Tdr = sylvop(ar, br, disc = true)
Tdrinv = invsylvop(ar,br, disc = true)
Tdrs = sylvop(as, bs, disc = true)
Tdrsinv = invsylvop(as,bs, disc = true)

Tdc = sylvop(ac, bc, disc = true)
Tdcinv = invsylvop(ac,bc, disc = true)
Tdcs = sylvop(acs, bcs, disc = true)
Tdcsinv = invsylvop(acs,bcs, disc = true)

@test check_ctranspose(Tdr) &&
      check_ctranspose(Tdrinv) &&
      check_ctranspose(Tdrs) &&
      check_ctranspose(Tdrsinv) &&
      check_ctranspose(Tdc) &&
      check_ctranspose(Tdcinv) &&
      check_ctranspose(Tdcs) &&
      check_ctranspose(Tdcsinv)

try
    T = invsylvop(as,bs,disc = true);
    T*rand(n*m);
    transpose(T)*rand(n*m);
    adjoint(T)*rand(n*m);
    T*rand(Float32,n*m);
    transpose(T)*rand(Float32,n*m);
    adjoint(T)*rand(Float32,n*m);
    T*ones(Int,n*m);
    transpose(T)*ones(Int,n*m);
    adjoint(T)*ones(Int,n*m);
    T = invsylvop(acs,bcs,disc = true);
    T*rand(n*m);
    T*complex(rand(n*m));
    transpose(T)*rand(n*m);
    transpose(T)*complex(rand(n*m));
    T'*rand(n*m);
    T'*complex(rand(n*m));
    @test true
catch
    @test false
end

@time x = Tdrinv*cr[:]
@test norm(Tdr*x[:]-cr[:])/norm(x[:]) < reltol

@time x = Tdrsinv*cr[:]
@test norm(Tdrs*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as,bs, disc = true)*cr[:]
@test norm(sylvop(as, bs, disc = true)*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as,bs', disc = true)*cr[:]
@test norm(sylvop(as, bs', disc = true)*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as',bs, disc = true)*cr[:]
@test norm(sylvop(as', bs, disc = true)*x[:]-cr[:])/norm(x[:]) < reltol

@time x = invsylvop(as',bs', disc = true)*cr[:]
@test norm(sylvop(as', bs', disc = true)*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as,bs, disc = true))*cr[:]
@test norm(transpose(sylvop(as, bs, disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as,bs', disc = true))*cr[:]
@test norm(transpose(sylvop(as, bs', disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as',bs, disc = true))*cr[:]
@test norm(transpose(sylvop(as', bs, disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as',bs', disc = true))*cr[:]
@test norm(transpose(sylvop(as', bs', disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as,bs, disc = true))*cr[:]
@test norm(adjoint(sylvop(as, bs, disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as,bs', disc = true))*cr[:]
@test norm(adjoint(sylvop(as, bs', disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as',bs, disc = true))*cr[:]
@test norm(adjoint(sylvop(as', bs, disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as',bs', disc = true))*cr[:]
@test norm(adjoint(sylvop(as', bs', disc = true))*x[:]-cr[:])/norm(x[:]) < reltol

@time x = Tdcinv*cc[:]
@test norm(Tdc*x[:]-cc[:])/norm(x[:]) < reltol

@time x = Tdcsinv*cc[:]
@test norm(Tdcs*x[:]-cc[:])/norm(x[:]) < reltol

@time x = Tdcsinv'*cr[:]
@test norm(Tdcs'*x[:]-cr[:])/norm(x[:]) < reltol

@time y = transpose(Tdrinv)*cr[:]
@test norm(transpose(Tdr)*y[:]-cr[:])/norm(y[:]) < reltol

@time y = transpose(Tdrsinv)*cr[:]
@test norm(transpose(Tdrs)*y[:]-cr[:])/norm(y[:]) < reltol

@time x = transpose(Tdrsinv)*cr[:]
@test norm(transpose(Tdrs)*x[:]-cr[:])/norm(x[:]) < reltol

@time y = Tdcinv'*cc[:]
@test norm(Tdc'*y[:]-cc[:])/norm(y[:]) < reltol

@time y = Tdcsinv'*cc[:]
@test norm(Tdcs'*y[:]-cc[:])/norm(y[:]) < reltol

x = rand(n*m);
@test abs(norm(Matrix(Tdr)'-Matrix(transpose(Tdr)))) < reltol
@test abs(norm(Matrix(Tdc)'-Matrix(Tdc'))) < reltol
@test norm(Matrix(Tdrinv)*Matrix(Tdr)-I) < reltol
@test norm(Matrix(Tdrsinv)*Matrix(Tdrs)-I) < reltol
@test norm(Matrix(Tdcinv)*Matrix(Tdc)-I) < reltol
@test norm(Matrix(Tdcsinv)*Matrix(Tdcs)-I) < reltol
@test norm(Tdrsinv'*x-invsylvop(as',bs', disc = true)*x) == 0
@test opnorm1(Tdrsinv-invsylvop(schur(ar),schur(br), disc = true)) == 0

@test sc*opnorm1(Tdr) < opnorm1est(Tdr)  &&
      sc*opnorm1(Tdrinv) < opnorm1est(Tdrinv)  &&
      sc*opnorm1(Tdrs) < opnorm1est(Tdrs)  &&
      sc*opnorm1(Tdrsinv) < opnorm1est(Tdrsinv)  &&
      sc*opnorm1(Tdc) < opnorm1est(Tdc)  &&
      sc*opnorm1(Tdcinv) < opnorm1est(Tdcinv)  &&
      sc*opnorm1(Tdcs) < opnorm1est(Tdcs)  &&
      sc*opnorm1(Tdcsinv) < opnorm1est(Tdcsinv)

@test sc*oprcondest(Tdr,Tdrinv) < 1/opnorm1(Tdr)/opnorm1(Tdrinv)  &&
      sc*oprcondest(Tdrs,Tdrsinv) < 1/opnorm1(Tdrs)/opnorm1(Tdrsinv)  &&
      sc*oprcondest(Tdc,Tdcinv) < 1/opnorm1(Tdc)/opnorm1(Tdcinv)  &&
      sc*oprcondest(Tdcs,Tdcsinv) < 1/opnorm1(Tdcs)/opnorm1(Tdcsinv) &&
      oprcondest(Tdr) == oprcondest(Tdr,Tdrinv) == oprcondest(Tdrinv)

@test opsepest(Tdrsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Tdrs)))  &&
      minimum(svdvals(Matrix(Tdrs))) <= sqrt(2)*n*opsepest(Tdrsinv)  &&
      opsepest(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.],disc=true)) == 0. &&
      opsepest(invsylvop([0. 1.; 0. 1.],-[0. 1.; 0. 1.],disc=true)) == 0.

end

# generalized Sylvester equations
@testset "Generalized Sylvester operators" begin

ar = rand(n,n)
br = rand(m,m)
cr = rand(n,n)
dr = rand(m,m)
er = rand(n,m)
ac = ar+im*rand(n,n)
bc = br+im*rand(m,m)
cc = cr+im*rand(n,n)
dc = dr+im*rand(m,m)
ec = er+im*rand(n,m)
as, cs = schur(ar,cr)
bs, ds = schur(br,dr)
acs, ccs = schur(ac,cc)
bcs, dcs = schur(bc,dc)

Tr = sylvop(ar, br, cr, dr)
Trinv = invsylvop(ar, br, cr, dr)
Trs = sylvop(as, bs, cs, ds)
Trsinv = invsylvop(as, bs, cs, ds)
Trs1 = sylvop(as, ds, cs, bs)
Trs1inv = invsylvop(as, ds, cs, bs)

Tc = sylvop(ac, bc, cc, dc)
Tcinv = invsylvop(ac, bc, cc, dc)
Tcs = sylvop(acs, bcs, ccs, dcs)
Tcsinv = invsylvop(acs, bcs, ccs, dcs)
@test size(Tr) == size(Tr') == size(Trinv')

@test check_ctranspose(Tr) &&
      check_ctranspose(Trinv) &&
      check_ctranspose(Trs) &&
      check_ctranspose(Trsinv) &&
      check_ctranspose(Tc) &&
      check_ctranspose(Tcinv) &&
      check_ctranspose(Tcs) &&
      check_ctranspose(Tcsinv)

try
    T = invsylvop(as,bs,cs,ds);
    T*rand(n*m);
    transpose(T)*rand(n*m);
    adjoint(T)*rand(n*m);
    T*rand(Float32,n*m);
    transpose(T)*rand(Float32,n*m);
    adjoint(T)*rand(Float32,n*m);
    T*ones(Int,n*m);
    transpose(T)*ones(Int,n*m);
    adjoint(T)*ones(Int,n*m);
    T = invsylvop(acs,bcs,ccs,dcs);
    T*rand(n*m);
    T*complex(rand(n*m));
    transpose(T)*rand(n*m);
    transpose(T)*complex(rand(n*m));
    T'*rand(n*m);
    T'*complex(rand(n*m));
    sylvop(ar,complex(br),cr,dr);
    invsylvop(ar,complex(br),cr,dr);
    @test true
 catch
    @test false
 end

@time x = Trinv*er[:];
@test norm(Tr*x-er[:])/norm(x[:]) < reltol

@time x = Trsinv*er[:]
@test norm(Trs*x[:]-er[:])/norm(x[:]) < reltol

@time x = Trs1inv*er[:]
@test norm(Trs1*x[:]-er[:])/norm(x[:]) < reltol

@time x = invsylvop(as, bs, cs, ds)*er[:]
@test norm(sylvop(as, bs, cs, ds)*x[:]-er[:])/norm(x[:]) < reltol

@time x = invsylvop(as', bs, cs', ds)*er[:]
@test norm(sylvop(as', bs, cs', ds)*x[:]-er[:])/norm(x[:]) < reltol

@time x = invsylvop(as, bs', cs, ds')*er[:]
@test norm(sylvop(as, bs', cs, ds')*x[:]-er[:])/norm(x[:]) < reltol

@time x = invsylvop(as', bs', cs', ds')*er[:]
@test norm(sylvop(as', bs', cs', ds')*x[:]-er[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as, bs, cs, ds))*er[:]
@test norm(transpose(sylvop(as, bs, cs, ds))*x[:]-er[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as, bs, cs, ds))*er[:]
@test norm(transpose(sylvop(as, bs, cs, ds))*x[:]-er[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as', bs, cs', ds))*er[:]
@test norm(transpose(sylvop(as', bs, cs', ds))*x[:]-er[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as, bs', cs, ds'))*er[:]
@test norm(transpose(sylvop(as, bs', cs, ds'))*x[:]-er[:])/norm(x[:]) < reltol

@time x = transpose(invsylvop(as', bs', cs', ds'))*er[:]
@test norm(transpose(sylvop(as', bs', cs', ds'))*x[:]-er[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as, bs, cs, ds))*er[:]
@test norm(adjoint(sylvop(as, bs, cs, ds))*x[:]-er[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as', bs, cs', ds))*er[:]
@test norm(adjoint(sylvop(as', bs, cs', ds))*x[:]-er[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as, bs', cs, ds'))*er[:]
@test norm(adjoint(sylvop(as, bs', cs, ds'))*x[:]-er[:])/norm(x[:]) < reltol

@time x = adjoint(invsylvop(as', bs', cs', ds'))*er[:]
@test norm(adjoint(sylvop(as', bs', cs', ds'))*x[:]-er[:])/norm(x[:]) < reltol

@time x = Tcinv*ec[:]
@test norm(Tc*x[:]-ec[:])/norm(x[:]) < reltol

@time x = Tcsinv*ec[:]
@test norm(Tcs*x[:]-ec[:])/norm(x[:]) < reltol

@time x = Tcsinv'*er[:]
@test norm(Tcs'*x[:]-er[:])/norm(x[:]) < reltol

@time y = transpose(Trinv)*er[:]
@test norm(transpose(Tr)*y[:]-er[:])/norm(y[:]) < reltol

@time y = transpose(Trsinv)*er[:]
@test norm(transpose(Trs)*y[:]-er[:])/norm(y[:]) < reltol

@time x = transpose(Trsinv)*er[:]
@test norm(transpose(Trs)*x[:]-er[:])/norm(x[:]) < reltol

@time x = Trs1inv'*er[:]
@test norm(Trs1'*x[:]-er[:])/norm(x[:]) < reltol

@time y = Tcinv'*ec[:]
@test norm(Tc'*y[:]-ec[:])/norm(y[:]) < reltol

@time y = Tcsinv'*ec[:]
@test norm(Tcs'*y[:]-ec[:])/norm(y[:]) < reltol

x = rand(n*m);
@test abs(norm(Matrix(Tr)'-Matrix(transpose(Tr)))) < reltol
@test abs(norm(Matrix(Tc)'-Matrix(Tc'))) < reltol
@test norm(Matrix(Trinv)*Matrix(Tr)-I) < reltol
@test norm(Matrix(Trsinv)*Matrix(Trs)-I) < reltol
@test norm(Matrix(Tcinv)*Matrix(Tc)-I) < reltol
@test norm(Matrix(Tcsinv)*Matrix(Tcs)-I) < reltol
@test norm(Trsinv'*x-invsylvop(as',bs', cs',ds')*x) == 0
@test opnorm1(Trsinv-invsylvop(schur(ar,cr),schur(br,dr))) == 0

@test sc*opnorm1(Tr) < opnorm1est(Tr)  &&
      sc*opnorm1(Trinv) < opnorm1est(Trinv)  &&
      sc*opnorm1(Trs) < opnorm1est(Trs)  &&
      sc*opnorm1(Trsinv) < opnorm1est(Trsinv)  &&
      sc*opnorm1(Tc) < opnorm1est(Tc)  &&
      sc*opnorm1(Tcinv) < opnorm1est(Tcinv)  &&
      sc*opnorm1(Tcs) < opnorm1est(Tcs)  &&
      sc*opnorm1(Tcsinv) < opnorm1est(Tcsinv)

@test sc*oprcondest(Tr,Trinv) < 1/opnorm1(Tr)/opnorm1(Trinv)  &&
      sc*oprcondest(Trs,Trsinv) < 1/opnorm1(Trs)/opnorm1(Trsinv)  &&
      sc*oprcondest(Tc,Tcinv) < 1/opnorm1(Tc)/opnorm1(Tcinv)  &&
      sc*oprcondest(Tcs,Tcsinv) < 1/opnorm1(Tcs)/opnorm1(Tcsinv) &&
      oprcondest(Tr) == oprcondest(Tr,Trinv) == oprcondest(Trinv)

@test opsepest(Trsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Trs)))  &&
      minimum(svdvals(Matrix(Trs))) <= sqrt(2)*n*opsepest(Trsinv)  &&
      opsepest(invsylvop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.])) == 0. &&
      opsepest(transpose(invsylvop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0. &&
      opsepest(adjoint(invsylvop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0. &&
      opsepest(invsylvop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.])) == 0. &&
      opsepest(transpose(invsylvop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0. &&
      opsepest(adjoint(invsylvop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0.

end

@testset "Sylvester system operators" begin

ar = rand(n,n)
br = rand(m,m)
cr = rand(n,n)
dr = rand(m,m)
er = rand(n,m)
fr = rand(n,m)
ac = ar+im*rand(n,n)
bc = br+im*rand(m,m)
cc = cr+im*rand(n,n)
dc = dr+im*rand(m,m)
ec = er+im*rand(n,m)
fc = fr+im*rand(n,m)
as, cs = schur(ar,cr)
bs, ds = schur(br,dr)
acs, ccs = schur(ac,cc)
bcs, dcs = schur(bc,dc)

Tr = sylvsysop(ar, br, cr, dr)
Trinv = invsylvsysop(ar, br, cr, dr)
Trs = sylvsysop(as, bs, cs, ds)
Trsinv = invsylvsysop(as, bs, cs, ds)

Tc = sylvsysop(ac, bc, cc, dc)
Tcinv = invsylvsysop(ac, bc, cc, dc)
Tcs = sylvsysop(acs, bcs, ccs, dcs)
Tcsinv = invsylvsysop(acs, bcs, ccs, dcs)
@test size(Tr) == size(Tr') == size(Trinv')

@test check_ctranspose(Tr) &&
      check_ctranspose(Trinv) &&
      check_ctranspose(Trs) &&
      check_ctranspose(Trsinv) &&
      check_ctranspose(Tc) &&
      check_ctranspose(Tcinv) &&
      check_ctranspose(Tcs) &&
      check_ctranspose(Tcsinv)

try
    T = invsylvsysop(as,bs,cs,ds);
    T*rand(2n*m);
    transpose(T)*rand(2n*m);
    adjoint(T)*rand(2n*m);
    T*rand(Float32,2n*m);
    transpose(T)*rand(Float32,2n*m);
    adjoint(T)*rand(Float32,2n*m);
    T*ones(Int,2n*m);
    transpose(T)*ones(Int,2n*m);
    adjoint(T)*ones(Int,2n*m);
    T = invsylvsysop(acs,bcs,ccs,dcs);
    T*rand(2n*m);
    T*complex(rand(2n*m));
    transpose(T)*rand(2n*m);
    transpose(T)*complex(rand(2n*m));
    T'*rand(2n*m);
    T'*complex(rand(2n*m));
    sylvsysop(ar,complex(br),cr,dr);
    invsylvsysop(ar,complex(br),cr,dr);
    @test true
 catch
    @test false
 end

@time xy = Trinv*[er[:];fr[:]]
@test norm(Tr*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@time xy = Trinv*[ec[:];fc[:]]
@test norm(Tr*xy-[ec[:];fc[:]])/norm(xy[:]) < reltol

@time xy = transpose(Trinv)*[er[:];fr[:]]
@test norm(transpose(Tr)*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@time xy = transpose(Trinv)*[ec[:];fc[:]]
@test norm(transpose(Tr)*xy-[ec[:];fc[:]])/norm(xy[:]) < reltol

@time xy = Tcinv*[ec[:];fc[:]]
@test norm(Tc*xy-[ec[:];fc[:]])/norm(xy[:]) < reltol

@time xy = Tcinv'*[ec[:];fc[:]]
@test norm(Tc'*xy-[ec[:];fc[:]])/norm(xy[:]) < reltol

@time xy = Trsinv*[er[:];fr[:]]
@test norm(Trs*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@time xy = Tcsinv*[er[:];fr[:]]
@test norm(Tcs*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@time xy = transpose(Trsinv)*[er[:];fr[:]]
@test norm(transpose(Trs)*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@time xy = Tcsinv*[er[:];fr[:]]
@test norm(Tcs*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@time xy = Tcsinv'*[er[:];fr[:]]
@test norm(Tcs'*xy-[er[:];fr[:]])/norm(xy[:]) < reltol

@test abs(norm(Matrix(Tr)'-Matrix(transpose(Tr)))) < reltol
@test abs(norm(Matrix(Tc)'-Matrix(Tc'))) < reltol
@test norm(Matrix(Trinv)*Matrix(Tr)-I) < reltol
@test norm(Matrix(Trsinv)*Matrix(Trs)-I) < reltol
@test norm(Matrix(Tcinv)*Matrix(Tc)-I) < reltol
@test norm(Matrix(Tcsinv)*Matrix(Tcs)-I) < reltol
@test opnorm1(Trsinv-invsylvsysop(schur(ar,cr),schur(br,dr))) == 0

@test sc*opnorm1(Tr) < opnorm1est(Tr)  &&
      sc*opnorm1(Trinv) < opnorm1est(Trinv)  &&
      sc*opnorm1(Trs) < opnorm1est(Trs)  &&
      sc*opnorm1(Trsinv) < opnorm1est(Trsinv)  &&
      sc*opnorm1(Tc) < opnorm1est(Tc)  &&
      sc*opnorm1(Tcinv) < opnorm1est(Tcinv)  &&
      sc*opnorm1(Tcs) < opnorm1est(Tcs)  &&
      sc*opnorm1(Tcsinv) < opnorm1est(Tcsinv)

@test sc*oprcondest(Tr,Trinv) < 1/opnorm1(Tr)/opnorm1(Trinv)  &&
      sc*oprcondest(Trs,Trsinv) < 1/opnorm1(Trs)/opnorm1(Trsinv)  &&
      sc*oprcondest(Tc,Tcinv) < 1/opnorm1(Tc)/opnorm1(Tcinv)  &&
      sc*oprcondest(Tcs,Tcsinv) < 1/opnorm1(Tcs)/opnorm1(Tcsinv) &&
      oprcondest(Tr) == oprcondest(Tr,Trinv) == oprcondest(Trinv)

@test opsepest(Trsinv)/n/sqrt(2) <= minimum(svdvals(Matrix(Trs)))  &&
      minimum(svdvals(Matrix(Trs))) <= sqrt(2)*n*opsepest(Trsinv)  &&
      opsepest(invsylvsysop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.])) == 0.   &&
      opsepest(transpose(invsylvsysop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0.  &&
      opsepest(adjoint(invsylvsysop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0.  &&
      opsepest(invsylvsysop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.])) == 0.  &&
      opsepest(transpose(invsylvsysop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0.  &&
      opsepest(adjoint(invsylvsysop([0. 1.; 0. 1.],[0. 1.; 0. 1.],-[0. 1.; 0. 1.],[0. 1.; 0. 1.]))) == 0.

end

end

end

max_stars_repo_name: andreasvarga/MatrixEquations

max_issues_repo_name: andreasvarga/MatrixEquations

max_forks_repo_name: andreasvarga/MatrixEquations


text:
# Visualization

function plot_state_city(state)
    qiskit.visualization.plot_state_city(state)
end

function plot_histogram(data)
    qiskit.visualization.plot_histogram(data)
end


text:
using Distributions

# Unit-tests for the VAR estimation

lags = 3
series = 6
length = 1000

#
# # Function for simulating Vector autoregressive process with given coefficient matrix
#

# This function simulates Vector autoregressive process taking those values as input:
# - length: length of the resulting time series
# - typ: type of the process, whether it should or should not have a constant, etc.
# - coefficients: matrix of coefficients with dimensions (series, series*lags) with ordering in each line as x_lag1, y_lag1, x_lag2, y_lag2, etc...
# - other_c: matrix of coefficients corresponding to potential constants and trends with size (series, ...) where ... is 0, 1, or 2, according to type

function VARsim(length, typ, coefficients, other_c)
	# Should be enough
	burnout = 200

	a, b = size(coefficients)
	lags = convert(Int, max(a,b)/min(a,b))
	series = min(a,b)

	typ=="None" && (other=fill(0.0, length + burnout))
	typ=="Const" && (other=fill(1.0, length + burnout))
	typ=="Trend" && (other=[1:(length+burnout)])
	typ=="Const and trend" && (other=[fill(1.0, length + burnout) [1:(length+burnout)]])

	y = fill(0.0, length + burnout, series)
	eps = rand(Normal(0, 1/100), (length + burnout, series))

	for i=(lags+1):(size(y)[1])
			y[i,:] = mapreduce((x) -> y[i-x,:]*coefficients[(1:series) + (x-1)*series, :], +, 1:lags) + eps[i,:] + other[i,:]*other_c'
	end

	y = y[(burnout+1):end,:]

	return y
end

#
# # Function simulating random VAR of given specification
#
# This function is suited for performing a unit-test of the VAR system.
# Arguments:
# - lags: number of lags in the system
# - series: number of series in the system
# - length: length of the resulting time series
# - typ: whether the system should not have constant, trend or have any combination of those.

function VARdatarandsim(lags, series, length, typ)
	typ=="None" && (other_c=fill(0.0, series))
	typ=="Const" && (other_c=rand(Uniform(), series))
	typ=="Trend" && (other_c=rand(Uniform(), series))
	typ=="Const and trend" && (other_c=[rand(Uniform(), series) rand(Uniform(), series)])

	# Set them so that it does not blow up!
	coefficients = rand(Uniform(0, 1/(lags*series)), (series, lags*series))'
	y = VARsim(length, typ, coefficients, other_c)
	return (y, [other_c', coefficients])
end

apply(vcat,[coeficients[(1:series) + (x-1)*series, :] for x=[lags:-1:1]]) - varEstimate(y, lags, "None").C

mapreduce((x)->VARdifcoef(lags, series, length, burnout), +, 1:60)/60


text:
using Base.Test
using Base.Threads
using Bio.Seq
using Bio.Intervals
using Libz
using DataFrames
using Gadfly

importall Bio.Intervals

include("pairs.jl")
include("energy.jl")
include("duplex.jl")
include("intervals.jl")
include("io.jl")
include("trie.jl")

function trie_profile( seq, range, num_steps; raw_output="" )
   stepsize = div(length(seq), num_steps)
   steps    = collect(50:stepsize:length(seq))
   build    = Vector{Float64}(length(steps))
   trav     = Vector{Float64}(length(steps))
   bulge    = Vector{Float64}(length(steps))
   mismatch = Vector{Float64}(length(steps))
   bulgemat = Vector{Float64}(length(steps))
   dep      = Vector{Int}(length(steps))
   cnt      = Vector{Int}(length(steps))
   int      = Vector{Int}(length(steps))

   df = DataFrame(length=[], buildtime=[], travtime=[], depth=[])

   if raw_output != ""
      const io = open( raw_output, "w" )
      const stream = ZlibDeflateOutputStream( io )
      const output = true
   else
      const output = false
   end

   for i in range
      for j in 1:length(steps)
         trie, build[j]   = @timed DuplexTrie{DNAAlphabet{2},UInt8}( seq[1:steps[j]], min(i, 5):i )
         cnt[j] = nodecount( trie.fwd )
         dep[j] = i
         res, trav[j]     = @timed traverse( trie, 1:10000 )
         int[j] = length(res)
         res, bulge[j]    = @timed traverse( trie, 1:10000, bulge_max=1 )
         res, mismatch[j] = @timed traverse( trie, 1:10000, mismatch_max=1 )
         res, bulgemat[j] = @timed traverse( trie, 1:10000, bulge_max=1, mismatch_max=1 )
      end
      if output
         for j in 1:length(steps)
            write( stream, string(i) * "\t" )
            write( stream, string(j) * "\t" )
            write( stream, string(steps[j]) * "\t" )
            write( stream, string(cnt[j]) * "\t" )
            write( stream, string(build[j]) * "\t" )
            write( stream, string(trav[j]) * "\t" )
            write( stream, string(bulge[j]) * "\t" )
            write( stream, string(mismatch[j]) * "\t" )
            write( stream, string(bulgemat[j]) * "\t" )
            write( stream, string(int[j]) * "\n" )
         end
      else
         df = join(df, DataFrame(length=steps, buildtime=build, travtime=trav, depth=dep), on=[:length, :buildtime, :travtime, :depth], kind = :outer)
      end
   end
   if !output
      buildplt = plot(df, x=:length, y=:buildtime, color=:depth, Geom.smooth)
      travplt  = plot(df, x=:length, y=:travtime,  color=:depth, Geom.smooth)

      return df,hstack(buildplt, travplt)
   else
      close(stream)
      close(io)
   end
end

function plot_profile( filename::String )
   df = readtable(filename, separator = '\t', header=false, names=[:depth, :step, :length, :nodecount, :buildtime,
                                                                   :travtime, :travbulgetime, :travmistime, :travbulgemistime, :duplexes])

   buildplt = plot(df, x=:length, y=:buildtime, color=:depth, Geom.smooth,
         Coord.cartesian(ymin=0),
         Guide.xlabel("Sequence Length"),
         Guide.ylabel("Trie Build Time (sec)"))
   nodeplt  = plot(df, x=:length, y=:nodecount,  color=:depth, Geom.line,
         Coord.cartesian(ymin=0),
         Guide.xlabel("Sequence Length"),
         Guide.ylabel("Node Count"))
   travplt  = plot(df, x=:length, y=:travtime,  color=:depth, Geom.smooth,
         Coord.cartesian(ymin=0),
         Guide.xlabel("Sequence Length"),
         Guide.ylabel("Traversal Time (sec)"))
  duplextravplt = plot(df, x=:duplexes, y=:travtime, color=:depth, Geom.line,
         Coord.cartesian(ymin=0),
         Guide.xlabel("Returned Duplexes"),
         Guide.ylabel("Traversal Time (sec)"))

   draw(PDF("profile.pdf", 9inch, 8inch), vstack(hstack(buildplt, nodeplt), hstack(travplt, duplextravplt)))

   res = stack(df, [:buildtime, :travtime, :travbulgetime, :travmistime, :travbulgemistime])
   timeplot = plot(res, x=:length, y=:value, color=:variable, Geom.smooth, Coord.cartesian(ymin=0, fixed=1.0),
                   Guide.xlabel("Sequence Length"),
                   Guide.ylabel("Traversal Time (sec)"))
   draw(PDF("profile_time.pdf", 4inch, 4inch), timeplot)
end


text:
const Point{T} = NTuple{2, T}

# ------------------------------------------------------------------------------------------
# Pedestrian
# ------------------------------------------------------------------------------------------
Base.@kwdef mutable struct Pedestrian{T<:Real}
    # mandatory properties
    id::Int
    pos::Point{T}
    vel::Point{T} = (0., 0.)
    target::Point{T} = pos

    # additional properties
    acc::Float64 = 0.5            # acceleration
    φ::Float64 = 3π/4             # maximum change of a pedestrian course
    radius::Float64 = 0.25        # social pedestrian size
    rlims::Point{T} = (0.1, 0.25) # physical pedestrian size

    # flags
    finished::Bool = false
    isexit::Bool = false

    # history
    position::Matrix{T} = zeros(0, 2)
    velocity::Matrix{T} = zeros(0, 2)
end

Pedestrian(id, pos; kwargs...) = Pedestrian(; id, pos, kwargs...)

"""
    isvalid(p::Pedestrian, pos, radius)

Checks wheter give position `pos` is valid with respect to pedestrian `p`.
"""
function isvalid(p::Pedestrian, pos::Point, radius::Real; social::Bool = false)
    radius_p = social ? p.rlims[2] : p.radius
    return distance(p.pos, pos) >= radius_p + radius
end

"""
    isvalid(p1::Pedestrian, p2::Pedestrian)

Checks wheter pedestrian `p2` has valid position with respect to pedestrian `p1`.
"""
function isvalid(p1::Pedestrian, p2::Pedestrian; kwargs...)
    return isvalid(p1, p2.pos, p2.radius; kwargs...)
end

# ------------------------------------------------------------------------------------------
# Room specification
# ------------------------------------------------------------------------------------------
abstract type RoomShape end
abstract type Obstacle end
abstract type Target end
abstract type Door <: Target end

struct Room
    shape::RoomShape
    obstacles::Vector{<:Obstacle}
    entrances::Vector{<:Door}
    exits::Vector{<:Door}
    checkpoints::Vector{<:Target}
end

"""
    isvalid(r::Room, pos::Point, r = 0)

Checks wheter give point `pos` is valid for given room configuration.
"""
function isvalid(room::Room, pos::Point, r = 0)
    # checks wheter given position lies inside the room (with zero radius)
    isvalid(room.shape, pos, 0) || return false

    # checks wheter given position is valid with respect to obstacles
    for o in room.obstacles
        isvalid(o, pos, r) || return false
    end

    # checks wheter given position lies inside the door area
    for d in vcat(room.entrances, room.exits)
        isvalid(d, pos, r) || return false
        isinside(d, pos, r) && return true
    end
    return isvalid(room.shape, pos, r)
end

# ------------------------------------------------------------------------------------------
# Model specification
# ------------------------------------------------------------------------------------------
abstract type TimeStrategy end
abstract type ScheduleStrategy end
abstract type TargetStrategy end
abstract type MoveStrategy end

Base.@kwdef struct Model{T<:Real}
    # room specification
    room::Room

    # strategies
    time_strategy::TimeStrategy
    schedule_strategy::ScheduleStrategy
    target_strategy::TargetStrategy
    move_strategy::MoveStrategy

    # pedestrians
    pedestrians::Dict{Int, Pedestrian} = Dict{Int, Pedestrian}()
    history::Dict{Int, Pedestrian} = Dict{Int, Pedestrian}()

    # additional parameters
    Δt::T = 0.05
    maxid::Ref{Int} = Ref(0)
    iter::Ref{Int} = Ref(0)
    safe_add::Bool = true
    stop::Ref{Bool} = Ref(false)
end

function Model(
    room,
    time_strategy,
    schedule_strategy,
    target_strategy,
    move_strategy,
    kwargs...
)
    return Model(;
        room,
        time_strategy,
        schedule_strategy,
        target_strategy,
        move_strategy,
        kwargs...
    )
end

function Base.show(io::IO, m::Model)
    println(io, "Pedestrian model:")
    println(io, " ⋅ room: ", m.room.shape)
    println(io, "   - checkpoints: ", length(m.room.checkpoints))
    println(io, "   - obstacles: ", length(m.room.obstacles))
    println(io, "   - entrances: ", length(m.room.entrances))
    println(io, "   - exits: ", length(m.room.exits))
    println(io, " ⋅ pedestrians: ", length(m.pedestrians))
    println(io, " ⋅ time strategy: ", typeof(m.time_strategy))
    println(io, " ⋅ schedule strategy: ", typeof(m.schedule_strategy))
    println(io, " ⋅ target strategy: ", typeof(m.target_strategy))
    println(io, " ⋅ move strategy: ", typeof(m.move_strategy))
    return
end

set_stop!(m::Model) = m.stop[] = true
maxid(m::Model) = m.maxid[]
npads(m::Model) = isempty(m.pedestrians)
iter(m::Model) = m.iter[]
allids(m::Model) = collect(keys(m.pedestrians))
allpedestrians(m::Model) = collect(values(m.pedestrians))

"""
    add_pedestrian!(m::Model, p::Pedestrian)

If the position of the pedestrian `p` is valid, adds `p` to the model `m` with smallest available `id`.
"""
function add_pedestrian!(m::Model, p::Pedestrian)
    if m.safe_add
        if !isvalid(m.room, p.pos, p.radius)
            @warn "invalid pedestrian with respect to room specification"
            return
        end
        for (~, p2) in m.pedestrians
            if !isvalid(p, p2)
                @warn "invalid pedestrian with respect to other pedestrians"
                return
            end
        end
    end

    m.maxid[] += 1
    p.id = m.maxid[]
    m.pedestrians[p.id] = p
    return
end

"""
    remove_pedestrian!(m::Model, p::Pedestrian)

Removes pedestrian `p` from the model `m`.
"""
function remove_pedestrian!(m::Model, p::Pedestrian)
    m.history[p.id] = deepcopy(p)
    delete!(m.pedestrians, p.id)
end

"""
    move_pedestrian!(m::Model, p::Pedestrian)

Move pedestrian `p` based on its position and velocity or remove `p` from model if in exit.
"""
function move_pedestrian!(m::Model, p::Pedestrian)
    p.velocity = vcat(p.velocity, [p.vel[1] p.vel[2]])
    p.position = vcat(p.position, [p.pos[1] p.pos[2]])
    p.finished && remove_pedestrian!(m, p)

    if p.isexit && norm(p.vel, 2)*m.Δt > distance(p.pos, p.target)
        p.finished = true
        remove_pedestrian!(m, p)
    else
        p.pos = p.pos .+ p.vel .* m.Δt
    end
end

function step!(m::Model)
    # add new pedestrians
    add_pedestrian!(m.time_strategy, m)

    # move all pedestrians
    for id in scheduler(m.schedule_strategy, m)
        pedestrian_step!(m.move_strategy, m, m.pedestrians[id])
    end

    # update counter
    m.iter[] += 1
    return
end

function run!(m::Model, iter)

    bar = Progress(iter; showspeed = true)
    for i in 1:iter
        step!(m)

        # progress bar
        k = length(keys(m.pedestrians))
        next!(bar; showvalues = [
                (:iter, i),
                (:pedestrians, k)
        ])
        if m.stop[]
            finish!(bar)
            break
        end
    end
    return
end