Implement a complex lattice to check scaling about lattices.jl

[Roger-luo]

@mbeach42 tried CardesianIndex for high dimensional lattice before #1 , I think we could just use this as site positions? Or do you think we need a type for site position ?

Update: Yeah, I think we need to consider implementing a unit cell interface or a lattice type called UnitCellLattice that construct a lattice from given unit cells, while its interface is still the same as others. I didn't need this feature myself when I implemented this package. This package was made for a Hamiltonian DSL here: https://github.com/Roger-luo/QuHamiltonian.jl the point is we can use lattice types as traits to dispatch other functions.

In order to allow extending the unit cells as custom type (since I suppose nobody want to write a long type parameter if the unit cell is too complicated), we need an extra type hierarchy for unit cells, e.g AbstractUnitCell. Then we will be able to dispatch different methods based on unit cells.

[carstenbauer]

What I mean there is we can always define an extra type hierarchy based on undirected lattice type by simply wrap it in a Directed struct to provide a default convention on the direction without much effort, since it is really cheap to have custom types in Julia, while having the same interface, just like what they do for graphs.

What I would like about this is that the user could step into the picture at whatever level of abstraction he wants: If he doesn't need labels, he gets a simpler lattice without labels.

The question to me is how much extra effort this would take. I reckon that it could be sometimes easier to immediately label sites during construction rather than taking an unlabeled lattice and adding labels retrospectively. But maybe should simply try and see?

[carstenbauer]

Another thing to note is that there are sometimes different unit cell definitions which produce the same lattice. Would be great if could support this, too (as LatticePhysics.jl) does, since this is relevant for a lot of simulations.

[carstenbauer]

Let me try to summarize the approaches taken (so far) by Lattices.jl and LatticePhysics.jl, they way I see it.

Lattices.jl:

• Every lattice is a type, e.g. Honeycomb{Parameters}, defined in a separate module, HoneycombLattice. Hence one can easily dispatch on a lattice.
• Different lattices should respect the AbstractLattice interface and thus should behave in the same way (duck typing). The interface is simple such that users could easily attach their own lattice implementations.
• There is no standard for the implementation of lattices.
• The available simple lattices are implemented as lazy lattices: iterators over sites and edges are calculated on-the-fly and there are no internal lists of all sites/edges.
• In the same way (and consequently), boundary conditions are imposed lazily as well (every lattice implements that on its own).
• Has no concept of a unit cell, Bravais lattice vector, etc.
• Lattices are immutable.

Status

• Available lattices: chain, square, (~ hyperrect)
• Userbase: presumably no one/very few
• Plotting: not available
• Pro: simple, rather idiomatic Julia
• Con: few lattices/functionality, very much early stage.

LatticePhysics.jl

• There are AbstractLattice, AbstractSite, etc. types. However, they are rather complex and interdependent. It's unlikely that other users will attach self-implemented concrete types.
• All lattices are implemented in the same way. Concrete Lattice, Unitcell, Site, etc. types are provided. One specifies a unit cell and constructs a lattice of a given size/shape/boundary condition out of it (see minimal example below).
• They all end up as parametric Lattice objects and the "which lattice am I" information is lost. We can't dispatch on different lattices but only their dimensionality etc.
• At its core, a lattice is just a list of sites and bonds (+ Bravais lattice vectors and unit cell information).
• Sites and edges, apart from the positional/directional information, have labels of arbitrary type, i.e. String, attached.
• Lattices are mutable, i.e. sites can be added/removed/relabeled etc. (see LatPhysLatticeModification.jl)
• Contains Luttinger-Tisza and Bandstructure calculation packages which are beyond the scope of a pure lattice library.

Minimal example

using LatPhysUnitcellLibrary, LatPhysLatticeConstruction
uc = getUnitcellHoneycomb()
lt = getLatticePeriodic(uc, 5)

# plotting
using Plots, LatPhysPlottingPlots
plot(lt)

Status

• Available lattices: (see LatPhysUnitcellLibrary.jl
  • 2D: square, kagome, square_octagon, triangular, honeycomb
  • 3D: cubic, diamond, fcc, hyperkagome, pryochlore, 8_3_a, 8_3_b, 8_3_c, 8_3_n, 9_3_a, 10_3_a, 10_3_b, 10_3_c, 10_3_d (see this paper for information on these numbered lattices).
• Userbase: ~ 6 users in our Cologne group (lattice definitions seem to be robust)
• Plotting:
  • LatPhysPlottingPyPlot
  • LatPhysPlottingPlots
  • SVG (not yet published but available offline)
• Pro: Lots of lattices
• Con:
  • complex (Users shouldn't care about unit cells etc. it should be a one-liner to get a lattice. Also, labels are irrelevant for most physics applications?)
  • not very idomatic Julia code (e.g. not duck typed)

Ok, this took a while, I'll comment on this later. Eventually we should try to get the best out of both worlds.

[ffreyer]

Also, labels are irrelevant for most physics applications?

The obvious use case would be the Kitaev model. Or any other model with edge-dependent interactions. It could also be (ab)used to cache previous results, which is what I'm doing in my classical Monte Carlo simulation for a 4-spin interaction. (Though I'm not using LatticePhysics for this)

[carstenbauer]

The obvious use case would be the Kitaev model. Or any other model with edge-dependent interactions.

Right. I've crossed out this line in my posting above.

The general question here is, what we mean by a lattice within the scope of this package.

• Merely an undirected graph that holds positional and connectivity information, or
• a directed graph with possibly extra information attached to sites and edges.

LatticePhysics.jl follows the second definition: for a 2x2 periodic square lattice it returns 16 directed edges (twice the number of undirected edges).

ALPS, in principle, also has a lattice as a directed graph with possible extra information attached to sites and edges (the xml definition files have edges with "source" and "target" information). However, if I create a 2x2 periodic square lattice I only get 8 directed edges (as for Square in Lattices.jl).

On a side note, while both LatticePhysics.jl and ALPS use linear indexes in bonds, Lattices.jl currently returns cartesian indexes:

julia> edges(Square(2,2)) |> collect
8-element Array{Tuple{Tuple{Int64,Int64},Tuple{Int64,Int64}},1}:
 ((1, 1), (2, 1))
 ((2, 1), (1, 1))
 ((1, 2), (2, 2))
 ((2, 2), (1, 2))
 ((1, 1), (1, 2))
 ((2, 1), (2, 2))
 ((1, 2), (1, 1))
 ((2, 2), (2, 1))

[Roger-luo]

As for directed graph, you can always get its edges by simply inverse the order of undirected graph

for (a, b) in edges(lattice)
     # edge ->
      a => b
     # edge <-
      b => a
end

The original design of this package is partly inspired by graph packages in JuliaGraphs. Since I think the way to implement lattices is indeed to treat them as special graphs with some annotation and symmetry, and we can always start with the simplest definition with undirect graph. I have a proposal for a refactor of this package based on @crstnbr and @ffreyer comments. To make a more idomatic Lattice package for Julians with potentially all feature requirements

On a side note, while both LatticePhysics.jl and ALPS use linear indexes in bonds, Lattices.jl currently returns cartesian indexes:

This is the default output of iterators, since that is what you got in the first place during calculation. But you can get its linear index by indexing the lattice with it, e.g lattice[site] gives you the id of this site.

Definition

Lattice here in this package defines a certain geometric object, which can be used to specify connections or interactions etc.

Type Hierarchies

• AbstractLattice{N} abstract type for all lattices

  • BoundedLattice{N, B} abstract lattices with boundaries, this is because we will also need things like random lattice (random connected interaction)
    • GraphLattice: a fallback type for general lattices which uses other graph package (AbstractGraph in https://github.com/JuliaGraphs) to store a general lattice
    • UnitCellLattice: a fallback type for general lattices (less general than GraphLattice) that assumes the lattice can be generated by given unit cell as subtype of AbstractUnitCell type. (it is basically the same Lattice type in LatticePhysics I suppose)
    • specialized lattice types (or just bindings to different UnitCellLattice with different UnitCell), e.g HyperRect (this should replace the Square and Chain in the implementation), Kagome, triangular etc.

• other lattices without geometric boundary

• AbstractUnitCell I haven't had a good idea what kind of type should this have yet, but potentially

  • UnitCell for the simplest version that only store the connections
  • MetaUnitCell which stores some extra info (e.g labels) in meta (similar to MetaGraph and SimpleGraph)

All the lattices type above assumes to represent a special kind of undirect graph, since we can always simply iterate through a geometric region (e.g edges, sites, neighbors) by iterate the same edge with different order twice. IIUC, I think the reason why ALPS or LatticePhysics uses
direct graph is mainly because you need to specify labels and bonds while constructing, which in original design you can't do it if the lattice is undirect (or the number of labels is not correct). But I think we should disentangle this two parts into different structure here. Labels is not the core functionality.

Why do we need specialized types? (Or why should we put some info in the type)

This is because these types can be used for further dispatching which enables people to write only a few lines with multiple dispatch. Moreover, in many cases, type information can also be used for high performance implementation. The implementation of Square and Chain are just an example of this idea, their performance is almost the same with manual straightforward version and have very little allocation during runtime.

But this is not useful if you want to iterate through one lattice over and over again. But I don't think this is an issue related to lattices itself,
it is a general issue for iterators, which means we can simply define a cached iterator (e.g cache(itr)), which will cache the collection of iterator as a Vector. We can let users to decide whether they want to cache or not explicitly. Thus this part is disentangled from the entire lattice implementation.

A practical reason for this is, if we store the bonds and sites (or just adjacency matrix) in our lattice type, when we calculate the faces or other region (in the future, e.g super edges as super graph), the iteration should still be cached, which make the implementation repeated. To follow the DRY (don't repeat yourself) principal, I think this should be separated.

FYI. Part of this infrastructure is implemented in Yao's CachedBlock as CacheServer.jl

Interfaces

• sites(lattice) : returns an iterable of all sites
• edges(lattice, distance=1): returns an iterable of all edges
• neighbors(lattice, site, distance=1): returns an iterable of neighbors of given site
• faces(lattice, distance=1) return an iterable of all faces

The interfaces mainly just ask for iteratbles, which make it very flexible and extensible.
We can return a cache(itr) by default as well.

[carstenbauer]

ping @janattig

[janattig]

First of all, thanks @crstnbr for making the nice comparison and bringing me into the discussion as well. Let me make some comments in the following as well:

The general question here is, what we mean by a lattice within the scope of this package.

• Merely an undirected graph that holds positional and connectivity information, or
• a directed graph with possibly extra information attached to sites and edges.

These are actually two different questions:

1. should the sites and bonds carry labels
2. should there be a convention of directionality

As for both my opinion is YES, since 1. you need to distinguish between different sites / bonds in nearly every physics case apart from lattices with only one site per unitcell and 2. most notably in the context of complex hopping / interaction as this introduces signs into the Hamiltonian matrix later on as is e.g. needed for Majorana models.

As for directed graph, you can always get its edges by simply inverse the order of undirected graph

for (a, b) in edges(lattice)
     # edge ->
      a => b
     # edge <-
      b => a
end

Although this is true, this is unfortunately not solving the problem of defining a convention of directionality.

[Roger-luo]

Thanks for the comment! this is very helpful, since I don't have as much experience on this direction as you. @ffreyer

However, I still have some question that I'm not sure, can you provide more description for convention of directionality? Do you mean
the lattice bonds/edges should have a default direction convention which can be used for later labelling?

should the sites and bonds carry labels

I didn't mean it shouldn't. But the abstraction is made at different level:

1. the most basic level, there is shouldn't be anything else except the geometric structure itself, this includes
   1.1 lattice plane, such as plane, sphere, etc.
   1.2 unit cells
2. the second level can be the labels

This design is actually used in many Julia packages already such as graphs. The LightGraphs.jl defines the most basic geometric structure of a graph, and then other packages will be built based on this package such as weighted graph (like you mentioned labels, weights are one kind of labels) SimpleWeightedGraphs.jl, StaticGraphs.jl, and if you like more generic labels, there's MetaGraph.

What I mean there is we can always define an extra type hierarchy based on undirected lattice type by simply wrap it in a Directed struct to provide a default convention on the direction without much effort, since it is really cheap to have custom types in Julia, while having the same interface, just like what they do for graphs.

Since lattices in principal is very similar to graphs (but not the same), I believe we can adopt similar design in the abstraction.

Moreover, this allows us to make use of Plot Receipts which is the official API of Plots instead of putting all the plot definition inside the core, which make it much easier to interact with different plotting backends.

[ffreyer]

I think you meant to ping @janattig.

[ffreyer]

Adding labes to sites or bonds can probably be done in a Node or Edge constructor respectively.

For example, we could have a NodeLattice{NodeType, EdgeType} calling NodeType(lazy_lattice, site_index, edge_list). The NodeType constructor should have enough information to generate labels and setup all required edges. We could also forward kwargs from the lattice constructor to the node/edge constructors to make it more generic.

Similarly, we can add an EdgeLattice calling the edge constructor first. I reckon we can use code generation for the different lattice geometries at this point. Or we could keep the simpler lattice type around and let people dispatch on that instead.

[janattig]

However, I still have some question that I'm not sure, can you provide more description for convention of directionality? Do you mean
the lattice bonds/edges should have a default direction convention which can be used for later labelling?

I mean that every edge i to j represents a matrix element H[i,j] so if these are complex, H[i,j] and H[j,i] are of opposite sign (since H is hermitian) and therefore one needs to specify which should carry which sign. This corresponds to fixing a directionality as in a bipartite lattice could be go from A to B always.

[janattig]

I think concerning the labeling it is much easier to provide sort of default labels when no label is needed by the user than to insert labels on a pre-existing structure later on.

[Roger-luo]

The question to me is how much extra effort this would take. I reckon that it could be sometimes easier to immediately label sites during construction rather than taking an unlabeled lattice and adding labels retrospectively. But maybe should simply try and see?

The extra effort is just one line @forward IIUC.

But I agree. I'm trying this in #17 , feel free to comment. I think the reason why doing this is because of the maintenance (flexibility for developers), as some of you might already notice, graph packages in JuliaGraphs are not developed by a single person, but a few people separately with a common interface, and this design will allow people to define their own staff separately, regarding to this package would eventually be some generic public package, it worth it to pay a bit more effort to make it easier to extend (flexible).

I think concerning the labeling it is much easier to provide sort of default labels when no label is needed by the user than to insert labels on a pre-existing structure later on.

Yes, this is exactly what EdgeLattice or whatever lattice type constructed on the basic geometric would do. The idea is inheritance over composition

struct BaseLattice <: AbstractLattice
end

struct LabelLattice <: AbstractLattice
     lattice::BaseLattice
     labels::LabelType
end

Although, actual implementation will be different from this, but this is just an equivalent way of implementing inheritance in Julia, this is equivalent to the following code in C++ (or Python)

class BaseLattice {

   # interface about geometry
}

class LabelLattice: public BaseLattice {
     LabelType labels

     # interface about the labels
}

The share the same interface for lattices, and in fact, nothing else except one line @forward is required to implement this. Thus to users
these two are basically looks the same since their behaviour is same except for the label interface.

Adding labes to sites or bonds can probably be done in a Node or Edge constructor respectively.

Yes, this is the cache iterator. In principal it should be equivalent to LatticePhysics's lattice in performance and storage complexity.