Commit a79b8a1b authored by Anton Akhmerov's avatar Anton Akhmerov
Browse files

implement consistent mode separation (closes #23)

parent aa96bfc0
Pipeline #4351 failed with stages
in 5 minutes and 47 seconds
......@@ -197,15 +197,8 @@ def order_schur(select, t, q, calc_ev=True, overwrite_tq=False):
trsen = getattr(lapack, ltype + "trsen")
# Figure out if select is a function or array.
isfun = isarray = True
try:
select(0)
except:
isfun = False
try:
select[0]
except:
isarray = False
isfun = callable(select)
isarray = hasattr(select, "__getitem__")
if not (isarray or isfun):
raise ValueError("select must be either a function or an array")
......
......@@ -14,6 +14,7 @@ from itertools import combinations_with_replacement
import numpy as np
import numpy.linalg as npl
import scipy.linalg as la
from scipy.spatial.distance import pdist, squareform
from .. import linalg as kla
from scipy.linalg import block_diag
from scipy.sparse import (identity as sp_identity, hstack as sp_hstack,
......@@ -393,7 +394,7 @@ def setup_linsys(h_cell, h_hop, tol=1e6, stabilization=None):
return Linsys(matrices, v_out, extract_wf)
def unified_eigenproblem(a, b=None, tol=1e6):
def unified_eigenproblem(a, b=None):
"""A helper routine for modes(), that wraps eigenproblems.
This routine wraps the regular and general eigenproblems that can arise
......@@ -406,60 +407,86 @@ def unified_eigenproblem(a, b=None, tol=1e6):
problem.
b : numpy array or None
The matrix on the right hand side of the generalized eigenvalue problem.
tol : float
The tolerance for separating eigenvalues with absolute value 1 from the
rest.
Returns
-------
ev : numpy array
An array of eigenvalues (can contain NaNs and Infs, but those
are not accessed in `modes()`) The number of eigenvalues equals
twice the number of nonzero singular values of
`h_hop` (so `2*h_cell.shape[0]` if `h_hop` is invertible).
evanselect : numpy integer array
Index array of right-decaying modes.
propselect : numpy integer array
Index array of propagating modes (both left and right).
vec_gen(select) : function
A function that computes the eigenvectors chosen by the array select.
ord_schur(select) : function
A function that computes the unitary matrix (corresponding to the right
eigenvector space) of the (general) Schur decomposition reordered such
that the eigenvalues chosen by the array select are in the top left
block.
An array of the eigenvalues of propagating modes.
The eigenvalues are sorted identically to ``prop_vecs``.
prop_vecs : numpy array
The array of eigenvectors corresponding to the propagating modes.
evan_vecs : numpy array
The array of Schur vectors corresponding to the right-decaying modes.
"""
if b is None:
eps = np.finfo(a.dtype).eps * tol
t, z, ev = kla.schur(a)
# Right-decaying modes.
select = abs(ev) > 1 + eps
# Propagating modes.
propselect = abs(abs(ev) - 1) < eps
vec_gen = lambda x: kla.evecs_from_schur(t, z, select=x)
ord_schur = lambda x: kla.order_schur(x, t, z, calc_ev=False)[1]
ev = 1 / ev.conj()
else:
eps = np.finfo(np.common_type(a, b)).eps * tol
s, t, z, alpha, beta = kla.gen_schur(a, b, calc_q=False)
# Right-decaying modes.
select = abs(alpha) > (1 + eps) * abs(beta)
# Propagating modes.
propselect = (abs(abs(alpha) - abs(beta)) < eps * abs(beta))
with np.errstate(divide='ignore', invalid='ignore'):
ev = alpha / beta
# Note: the division is OK here, since we later only access
ev = (beta / alpha).conj()
# Note: the division is OK here, since we later only access the
# eigenvalues close to the unit circle
vec_gen = lambda x: kla.evecs_from_gen_schur(s, t, z=z, select=x)
ord_schur = lambda x: kla.order_gen_schur(x, s, t, z=z,
calc_ev=False)[2]
return ev, select, propselect, vec_gen, ord_schur
# The eigenvalues either form pairs (lambda, 1/lambda^*) or they are lying
# on the unit circle. Here we form a good quality splitting of modes into
# pairs of approximately inverses and single modes close to the unit
# circle. The most tricky case is when we have a group of propagating modes
# with the same momentum (due to a level crossing or even a conservation
# law). The errors would randomly move these modes around a bit, and can
# result in the two eigenvalues being closer to each other's inverses than
# the unit circle. Since estimating the error is hard (we're solving a
# potentially badly conditioned generalized eigenvalue problem), we need to
# use a criterion for dealing with this case that isn't tied to a specific
# precision.
# We need to figure out when an eigenvalue computed with an error is more
# likely to be its own partner than a partner of a nearby eigenvalue. To
# avoid false negatives we declare modes propagating if there is no partner
# closer than the *square* of the distance from the unit circle.
outer = abs(ev) > 1
inner = np.logical_not(outer)
ev[outer] = 1./ev[outer].conj()
distances_unit_circle = 1 - abs(ev)
coords = np.vstack([ev.real, ev.imag])
distances = squareform(pdist(coords.T))
distances[np.ix_(inner, inner)] = distances[np.ix_(outer, outer)] = np.inf
distances[np.diag_indices_from(distances)] = distances_unit_circle**2
partners = np.empty((len(ev),))
exact_matches = np.where(distances_unit_circle == 0)[0]
partners[exact_matches] = exact_matches
paired = set(exact_matches)
for index in np.argsort(distances_unit_circle)[::-1]:
if index in paired:
continue
partner = np.argmin(distances[index])
paired.add(partner)
distances[:, partner] = np.inf
distances[:, index] = np.inf
partners[index] = partner
partners[partner] = index
# Propagating modes.
propagating = (partners == np.arange(len(ev)))
# Right-decaying modes.
decaying = inner & ~propagating
print(ev - np.sign(ev))
print(inner)
print(decaying)
print((ev - np.sign(ev))[decaying])
return (ev[propagating], vec_gen(propagating),
ord_schur(decaying)[:, :np.sum(decaying)])
def phs_symmetrization(wfs, particle_hole):
......@@ -770,6 +797,7 @@ def make_proper_modes(lmbdainv, psi, extract, tol, particle_hole,
full_psi[:, indx[len(indx)//2:]] = out
psi[:, indx[len(indx)//2:]] = psi[:, indx[len(indx)//2:]].dot(rot)
print(velocities)
if np.any(abs(velocities) < vel_eps):
raise RuntimeError("Found a mode with zero or close to zero velocity.")
if 2 * np.sum(velocities < 0) != len(velocities):
......@@ -853,29 +881,19 @@ def compute_block_modes(h_cell, h_hop, tol, stabilization,
# Defer most of the calculation to helper routines.
matrices, v, extract = setup_linsys(h_cell, h_hop, tol, stabilization)
ev, evanselect, propselect, vec_gen, ord_schur = unified_eigenproblem(
*(matrices + (tol,)))
ev, prop_vecs, evan_vecs = unified_eigenproblem(*(matrices))
# v is never None.
# h_hop.shape[0] and v.shape[1] not always the same.
# Adding this makes the difference for some tests
n = v.shape[1]
nrightmovers = np.sum(propselect) // 2
nevan = n - nrightmovers
evan_vecs = ord_schur(evanselect)[:, :nevan]
# Compute the propagating eigenvectors.
prop_vecs = vec_gen(propselect)
# Compute their velocity, and, if necessary, rotate them.
# prop_vecs here is 'psi' in make_proper_modes, i.e. the wf in the SVD
# basis. It is in turn used to construct vecs and vecslmbdainv (the
# propagating parts). The evanescent parts of vecs and vecslmbdainv come
# from evan_vecs above.
prop_vecs, real_space_data = make_proper_modes(ev[propselect], prop_vecs,
extract, tol, particle_hole,
time_reversal, chiral)
# Transform the propagating eigenstates to be eigenstates of current with
# eigenvalues ±1.
prop_vecs, real_space_data = make_proper_modes(ev, prop_vecs, extract, tol,
particle_hole=particle_hole,
time_reversal=time_reversal,
chiral=chiral)
vecs = np.c_[prop_vecs[n:], evan_vecs[n:]]
vecslmbdainv = np.c_[prop_vecs[:n], evan_vecs[:n]]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment