RFC: MeshSpec and linearmesh, with support for bandcuts (i.e. lifting)

[pablosanjose]

Closes: #64

This implements the proposal by @BacAmorim in #64 that unifies marchingmesh and linearmesh to both construct subtypes of MeshSpec. This kind of type holds information necessary to build a Mesh and a lift function, once h is provided. Multiple dispatch on functions buildmesh(spec, br) and buildlift(spec, br) take care of that, where br is the Bravais matrix of h (which is actually the only thing required to build the mesh from the spec).

We thus now have a new API for bandstructures, which is the one a typical user is likely to use, based on MeshSpec objects

bandstructure(::Union{Hamiltonian,ParametricHamiltonian},  ::MeshSpec; kw...)

where we currently have two types of MeshSpec, MarchingMeshSpec and LinearMeshSpec, constructed using marchingmesh(ranges; resolution, axes) and linearmesh(nodes; resolution, samelength, closed) functions.

The linearmesh function implements named points in the Brillouin zone as nodes for polygonal cuts. A Γ-K-M-Γ cut, for example, would be done like

bandstructure(h, linearmesh(:Γ, :K, :M, :Γ; resolution = (10, 6, 10)))

or alternatively specifying one or more of the points as phi tuples,

bandstructure(h, linearmesh((0,0), (2pi/3, 2pi/3), (2pi/3, 0), (0,0); resolution = (10, 6, 10)))

Note that we can specify different resolutions per segment in the cut, or an equal resolution with a single integer.

The updated marchingmesh::MeshSpec API has a syntax based on tuples, instead of ranges, in keeping with the above

bandstructure(h, marchingmesh((-pi, pi), (0, 2pi); resolution = (10, 20)))

The old API bandmesh(h, mesh; lift = ..., kw...) is still available. The MeshSpec API above actually calls bandmesh(h, mesh; lift = ..., kw...) under the hood after computing the mesh and the lift.

[pablosanjose]

The last commit includes a fix in resolve_degeneracies!, which before would only explore finite-differences along mesh directions. The most robust way is to use finite difference along Brillouin zone directions, even if they lie outside of the mesh. In this way we fix cases which would be perfectly degenerate along a cut, but not in orthogonal directions.

This however has broken some tests with ParametricHamiltonian. Looking into ways to fix that. [EDIT: fixed, hopefully]

[codecov-commenter]

Codecov Report

Merging #70 into master will increase coverage by 4.65%.
The diff coverage is 92.06%.

@@            Coverage Diff             @@
##           master      #70      +/-   ##
==========================================
+ Coverage   56.23%   60.88%   +4.65%     
==========================================
  Files          15       15              
  Lines        2287     2572     +285     
==========================================
+ Hits         1286     1566     +280     
- Misses       1001     1006       +5     

Impacted Files	Coverage Δ
src/Quantica.jl	100.00% <ø> (ø)
src/plot_vegalite.jl	0.00% <0.00%> (ø)
src/diagonalizer.jl	56.90% <88.46%> (+4.06%)	⬆️
src/mesh.jl	93.91% <93.42%> (+0.32%)	⬆️
src/bandstructure.jl	92.30% <94.44%> (-0.32%)	⬇️
src/parametric.jl	82.75% <100.00%> (+1.14%)	⬆️
src/tools.jl	63.02% <100.00%> (+1.31%)	⬆️
src/hamiltonian.jl	71.20% <0.00%> (+3.63%)	⬆️
... and 1 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1f85af0...cea2477. Read the comment docs.

[pablosanjose]

Updated the OP with more details.

[BacAmorim]

I think the general API bandstructure(::Hamiltonian, ::MeshSpec) makes sense for the me (I do not know if Me <: AbstractUser). I just have a few comments:

• the keyword resolution is a bit misleading. It reminds of specifying the tolerance in an adaptive integration routine (a bit like PrecisionGoal and AccuracyGoal in Mathematica), when in fact it is just the number of points. I think a keyword npoints would be more descriptive.
• If I understood correctly the code, if in linearmesh we set resolution = 10, this means that each segment of the linearmesh will have 10 points. I think it would be more natural for resolution to specify the total number of points in the path instead. (I am aware that in marchingmesh, resolution = 10 would mean that the mesh has 10 points along each dimension)
• In my mind it stills makes more sense to include lift has an additional field in MarchingMeshSpec. I think that the meshspec should contain the information of the space we are actually sampling, such that buildlift(::MeshSpec) should always return the lift function that is being used. I think it is more natural to write

b = bandstructure(h, marchingmesh((0, 2π), npoints = 100, lift = φ -> (φ, -φ)))

• Why can't we use a lift function with bandstructure(matrixf::Function, mesh::Mesh; kw...)? I think it makes as much sense as for an Hamiltonian
• In the documentation of bandstructure(matrixf::Function, mesh::Mesh; kw...), we can read "Note that ϕ is either a Tuple or an SVector,...". Wouldn't it be better to have a unique form of internally storing ϕ?

[pablosanjose]

the keyword resolution is a bit misleading. It reminds of specifying the tolerance in an adaptive integration routine (a bit like PrecisionGoal and AccuracyGoal in Mathematica), when in fact it is just the number of points. I think a keyword npoints would be more descriptive.

I think I agree. Let's settle on points.

If I understood correctly the code, if in linearmesh we set resolution = 10, this means that each segment of the linearmesh will have 10 points. I think it would be more natural for resolution to specify the total number of points in the path instead. (I am aware that in marchingmesh, resolution = 10 would mean that the mesh has 10 points along each dimension)

This is problematic, because we want the provided nodes to be included in the mesh. If resolution - 1 is not a multiple of the number of segments that will not be the case.

In my mind it stills makes more sense to include lift has an additional field in MarchingMeshSpec. I think that the meshspec should contain the information of the space we are actually sampling, such that buildlift(::MeshSpec) should always return the lift function that is being used. I think it is more natural to write

This is again problematic. The lift function depends on the dimensions of the hamiltonian, so it cannot be an intrinsic part of the mesh. It is literally the mapping between the mesh and the Hamiltonian parameter space, so it needs to go outside.

Why can't we use a lift function with bandstructure(matrixf::Function, mesh::Mesh; kw...)? I think it makes as much sense as for an Hamiltonian

Because in Julia we cannot know the number of arguments matrixf accepts. That is a consequence of multiple dispatch: a function can have several methods, each with a different signature

n the documentation of bandstructure(matrixf::Function, mesh::Mesh; kw...), we can read "Note that ϕ is either a Tuple or an SVector,...". Wouldn't it be better to have a unique form of internally storing ϕ?

Thanks for looking into the docs! Internally almost everything (including the vertices of a mesh) are StaticArrays. It is just that, for convenience, we accept tuples that we convert later (they are much easier to type)

[BacAmorim]

This is problematic, because we want the provided nodes to be included in the mesh. If resolution-1 is not a multiple of the number of segments that will not be the case.

That is true. But we can work around that by defining a function that distributes the number of points has equally as possible. Something like:

function points_per_segment(s::LinearMeshSpec{N,L,T,R}, br) where {N, L, T, R<:Integer}
    
	nverts = length(s.vertices)
	npts = only(s.resolution)
	seg_lengths = segment_lengths(s, br) # are nomalized to sum(seg_lengths) == 1
	seg_npts =  [floor(Int, (npts-1)*seg_lengths[n]) for n in 1:(nverts-1)] # underestimate number of points per segment
	seg_remainder = [npts*seg_lengths[n] - seg_npts[n] for n in 1:(nverts-1)] # how many points can we still fit in each segment
	perm = sortperm(seg_remainder, rev = true) # sort segments by the number of points we can still fit in them
	r = npts - sum(seg_npts) # number of points that haven't been distributed yet
	p = 1
	while r>0
		seg_npts[perm[p]] += 1
		r -= 1
		p += 1
	end
	
	return ntuple(i->seg_npts[i], Val(N-1))
    
end   

function points_per_segment(s::LinearMeshSpec{N,L,T,R}, br) where {N, L, T, R} # all other cases
    
    return s.resolution
    
end  

Then we would just need to modify:

function buildmesh(s::LinearMeshSpec{N}, br) where {N}
    ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), points_per_segment(s, br)) # modify this line
    verts = SVector.(first(ranges))
    for r in Base.tail(ranges)
        pop!(verts)
        append!(verts, SVector.(r))
    end
    s.closed && pop!(verts)
    nv = length(verts)
    nodefunc = idx_to_node(s, br)
    verts .= nodefunc.(verts)
    adjmat = sparse(vcat(1:nv-1, 2:nv), vcat(2:nv, 1:nv-1), true, nv, nv)
    s.closed && (adjmat[end, 1] = adjmat[1, end] = true)
    return Mesh(verts, adjmat)
end

and I think everything should work.

This is again problematic. The lift function depends on the dimensions of the hamiltonian, so it cannot be an intrinsic part of the mesh. It is literally the mapping between the mesh and the Hamiltonian parameter space, so it needs to go outside.

The lift function has to be compatible with the dimensions of the Hamiltonian, yes. But I don't understand how it 'depends' on the dimensions of the Hamiltonian. What would be the difference in writing (example in the docs):

bandstructure(h, marchingmesh((0, 2π); points = 25); lift = φ -> (φ, 0))

or

bandstructure(h, marchingmesh((0, 2π); points = 25, lift = φ -> (φ, 0)))

with lift stored in a field MarchingMeshSpec.lift = lift, such that buildlift(s, br) = s.lift? In both cases the user is responsible in making sure that the lift function is compatible with the dimensions of the Hamiltonian and the mesh. For the same reason, I don't understand the next point:

Because in Julia we cannot know the number of arguments matrixf accepts. That is a consequence of multiple dispatch: a function can have several methods, each with a different signature

I understand this, but I do not see how the picture is different if we just provide the mesh, instead of mesh+lift. How do we know if each vertex in the mesh provided by a user has the correct number of elements that are taken as an argument by matrixf? We do not know. It us up to the user to make sure that the mesh he provides is compatible with the function matrixf. In the same maner, if there were a lift function, it would still be up to the user to make sure that the output of the lift function is compatible with matrixf (and with the mesh vertices). In both cases, it is up to the user to make sure that things are compatible. I don't see how we can protect him from this. If he provides an incompatible mesh/lift we will get an error saying that there is no method of matrixf compatible.

Now these two points are not very important. You have already made the functionality available. I was just trying to understand the rational for the way you implemented it. (I always learn new Julia things everytime I look at your code!)

EDIT: fiex typos in points_per_segment

[BacAmorim]

I found a possible problem. If I try:

a1 = SVector(0.5, sqrt(3)/2)
a2 = SVector(-0.5, sqrt(3)/2)

lat = lattice(
    bravais(a1, a2), 
    sublat(SVector(0.0, 0.0), name = :A), 
    sublat(SVector(0.0, 1/sqrt(3)), name = :B)
)
graphene = hamiltonian(
    lat, 
    hopping(-1.0, sublats = :A => :B, range = 1) + 
    hopping(-1.0, sublats = :B => :A, range = 1)
)

linearspec = linearmesh(:Γ, :K, :M, :Γ; resolution = (5, 2, 8)) 

I would expect the mesh generated from the spec to have 5 + 2 + 8 = 15 points
But I actually get:

buildmesh(linearspec, bravais(graphene)

Mesh{1}: mesh of a 1-dimensional manifold
  Vertices   : 13
  Edges      : 12

[BacAmorim]

If we want the case with resolution = number to specify the total number of points in the path, we can define a function:

Suggested change

	function points_per_segment(s::LinearMeshSpec{N,L,T,R}, br) where {N, L, T, R<:Integer}
	nverts = length(s.vertices)
	npts = only(s.resolution)
	seg_lengths = segment_lengths(s, br) # lengths are nomalized to sum(seg_lengths) == 1
	seg_npts = [floor(Int, (npts-1)*seg_lengths[n]) for n in 1:(nverts-1)] # underestimate number of points per segment
	seg_remainder = [npts*seg_lengths[n] - seg_npts[n] for n in 1:(nverts-1)] # how many points can we still fit in each segment
	perm = sortperm(seg_remainder, rev = true) # sort segments by the number of points we can still fit in them
	r = npts - 1 - sum(seg_npts) # number of points that haven't been distributed yet
	p = 1
	while r>0
	seg_npts[perm[p]] += 1
	r -= 1
	p += 1
	end
	# for convinince in buildmesh: sum(npts) == number of desired points - 1. Ensures the correct numbed of points
	return ntuple(i->seg_npts[i], Val(N-1))
	end
	function points_per_segment(s::Quantica.LinearMeshSpec{N,L,T,R}, br) where {N, L, T, R}
	npts = [n for n in s.resolution]
	npts[end] -= 1 # for convinince in buildmesh. sum(npts) == number of desired points - 1. Ensures the correct numbed of points
	return ntuple(i->npts[i], Val(N-1))
	end

[BacAmorim]

With the previous definition of points_per_segment, buildmesh can be written as (such that the number of vertices in the returned mesh ==sum(resolution)):

Suggested change

	function buildmesh(s::LinearMeshSpec{N}, br) where {N}
	ranges = ((i, r) -> range(i, i+1, length = r)).(ntuple(identity, Val(N-1)), s.points)
	verts = SVector.(first(ranges))
	for r in Base.tail(ranges)
	pop!(verts)
	append!(verts, SVector.(r))
	end
	s.closed && pop!(verts)
	nv = length(verts)
	nodefunc = idx_to_node(s, br)
	verts .= nodefunc.(verts)
	adjmat = sparse(vcat(1:nv-1, 2:nv), vcat(2:nv, 1:nv-1), true, nv, nv)
	s.closed && (adjmat[end, 1] = adjmat[1, end] = true)
	return Mesh(verts, adjmat)
	end
	function buildmesh(s::Quantica.LinearMeshSpec{N}, br) where {N}
	nsegs = length(s.vertices)-1
	seg_npts = points_per_segment(s, br)
	nodes =linearmesh_nodes(s, br)
	verts = [SVector(nodes[1])]
	for n in 1:nsegs
	for i in 1:seg_npts[n]
	push!(verts, SVector(nodes[n] + (nodes[n+1]-nodes[n])*i/seg_npts[n]))
	end
	end
	nv = length(verts)
	adjmat = sparse(vcat(1:nv-1, 2:nv), vcat(2:nv, 1:nv-1), true, nv, nv)
	s.closed && (adjmat[end, 1] = adjmat[1, end] = true)
	return Mesh(verts, adjmat)
	end

[pablosanjose]

That is true. But we can work around that by defining a function that distributes the number of points has equally as possible.

I'm afraid I see that as trying to be too clever for our own good :-D. I'd rather opt for behaviour that is more predictable and easier to maintain and understand

The lift function has to be compatible with the dimensions of the Hamiltonian, yes. But I don't understand how it 'depends' on the dimensions of the Hamiltonian. What would be the difference in writing (example in the docs):

It is possible, but it is not conceptually cleaner. Again, why would you want to mix the concepts of a mesh with the mapping of that mesh to a Hamiltonian? In fact, if you do that, how can a LinearMeshSpec compute its appropriate lift function for a given h without knowing the bravais(h). Note that with this desing, something like linearmesh(pts...; points = 25, lift = φ -> (φ, 0)) should be able to instantiate a MarchingMeshSpec object, independently of h.

How do we know if each vertex in the mesh provided by a user has the correct number of elements that are taken as an argument by matrixf?

You've got a point. I need to think about this again, as I feel there was a more powerful reason for this decision (it's been a while)

I would expect the mesh generated from the spec to have 5 + 2 + 8 = 15 points

Note that 5 is the points of a segment including endpoints (see docs). Hence, it is actually 4+1+8 = 13. The code is working as intended.

[BacAmorim]

I'm afraid I see that as trying to be too clever for our own good :-D. I'd rather opt for behaviour that is more predictable and easier to maintain and understand

Fair enough!

Again, why would you want to mix the concepts of a mesh with the mapping of that mesh to a Hamiltonian?

It is just that right now, for a LinearMeshSpec the lift function is built from it (and, yes, from the Bravais basis) using buildlift. For me it seems natural to have a MeshSpec to do two things: (i) generate a mesh in some space, (ii) provide a way to map the mesh to the target space of Hamiltonian arguments (by possibly supplying the bravais basis of the Hamiltonian). This is what LinearMeshSpec does right now. However, it is not what MarchingMeshSpec does, for which any (non-trivial) lift functions is provided outside the spec. It seems two ways of doing similar things.

Note that 5 is the points of a segment including endpoints (see docs). Hence, it is actually 4+1+8 = 13. The code is working as intended.

Oh sorry, I missed that in the docs!

[pablosanjose]

It is just that right now, for a LinearMeshSpec the lift function is built from it (and, yes, from the Bravais basis) using buildlift. For me it seems natural to have a MeshSpec to do two things: (i) generate a mesh in some space, (ii) provide a way to map the mesh to the target space of Hamiltonian arguments (by possibly supplying the bravais basis of the Hamiltonian). This is what LinearMeshSpec does right now. However, it is not what MarchingMeshSpec does, for which any (non-trivial) lift functions is provided outside the spec. It seems two ways of doing similar things.

But it is not different, really. All MeshSpecs have a default lift (generated by buildlift(spec, bravais)) that is used if no lift is provided by the user. In the case of a LinearMeshSpec it is the polygonal function across the nodes. In the case of a MarchingMeshSpec it is the identity (if the dimension of the mesh matches the Hamiltonian's), a padding with zeros of each vertex (if the dimension is lower) or an projection-out of the excess dimensions (if the dimension is higher). You may argue that these are trivial lifts, and they are, but it is what makes sense. Note again that in both the MarchingMeshSpec and LinearMeshSpec cases, a different lift function can be provided by the user.

How do we know if each vertex in the mesh provided by a user has the correct number of elements that are taken as an argument by matrixf?

You are right about lifts and matrixf::Function, there is no real reason to disallow them! I'm pushing a commit with this generalization in a bit.

[pablosanjose]

What I think we must definitely disallow is the method bandstructure(::Function, ::MeshSpec;...), because in general we need a Bravais matrix to build the mesh and the lift. Trying to make it work in particular cases where it may make sense is a bad idea.

[pablosanjose]

One more stylistic tweak: since buildmesh(s, br) and buildlift(s, br) are part of the user-facing exported API, they should not require a br::SMatrix Bravais matrix, since the fact that we use SMatrix for Bravais matrices is an implementation detail. To avoid leaking it out, we now document only buildmesh(::MeshSpec, ::Union{Hamiltonian,ParametricHamiltonian}), and similarly for buildlift. This should be the methods used by a user. Another reason for this is that a different, future ::MeshSpec type that we add could require more info about h than just its Bravais matrix.

I plan to merge this soon. Any more suggestions?

[BacAmorim]

Note again that in both the MarchingMeshSpec and LinearMeshSpec cases, a different lift function can be provided by the user.

Ok, this convinced me. So the idea is that each MeshSpec has one default lift function. The user can overwrite this, but that is outside the spec. I see the logic now.

Regarding bandstructure(::Function, ::Mesh;...), well the lift can be simply implemented by the user by making

bandstructure(x->matrixf(lift(x)), mesh)

So maybe that was your initial reasoning?

[pablosanjose]

Ok, this convinced. So the idea is that each MeshSpec has one default lift function. The user can overwrite this, but that is outside the spec. I see the logic now.

To be precise, one MeshSpec has one default lift function for a given Hamiltonian, but yes. Note that a marchingmesh can be used with a non-default lift function (and this is likely to be actually useful), for example to make a 2D cut along an arbitrary plane in a 3D bandstructure. Here, then, marchingmeh is the 2D plane, and the lift dictates how we cut the 3D Brillouin zone with it.

So maybe that was your initial reasoning?

Well, that is how it works now. Any bandstructure(::Hamiltonian...) call ends up in a bandstructure(::Function...) call, which encodes the lift precisely like that. But what you made me realize is that that is not a reason to exclude a lift option when the user calls bandstructure(f::Function...) directly, which she could want to do in cases where the bandstructure of a non-Hamiltonian object is desired (this was part of the generalization in #17). In this case, this PR merely wraps the lift into a new function f´ = phi -> f(applylift(f, phi)), and calls _bandstructure(f´,...) instead, which is as simple as it can get.

EDIT: I had misunderstood your comment. Yes, precisely, that's what's done in commit 4aa3647