From 64d9d4be7e193481352987bac1905b9bfe313097 Mon Sep 17 00:00:00 2001 From: Joseph Weston <joseph@weston.cloud> Date: Mon, 9 Sep 2019 13:35:27 +0200 Subject: [PATCH] update abstract System to support vectorization Add the 'subgraphs', 'terms', and 'site_arrays' attributes, and 'hamiltonian_term' method and implement 'hamiltonian_submatrix', 'cell_hamiltonian' and 'inter_cell_hopping' in a vectorized way. --- doc/source/reference/kwant.system.rst | 13 +- kwant/_system.pyx | 386 +++++++++++++++++++++++++- kwant/system.py | 166 +++++++++-- 3 files changed, 534 insertions(+), 31 deletions(-) diff --git a/doc/source/reference/kwant.system.rst b/doc/source/reference/kwant.system.rst index 357765c9..2a9bd89c 100644 --- a/doc/source/reference/kwant.system.rst +++ b/doc/source/reference/kwant.system.rst @@ -29,11 +29,22 @@ Sites SiteFamily Systems --------- +------- .. autosummary:: :toctree: generated/ System + VectorizedSystem InfiniteSystem + InfiniteVectorizedSystem FiniteSystem + FiniteVectorizedSystem PrecalculatedLead + +Mixin Classes +------------- +.. autosummary:: + :toctree: generated/ + + FiniteSystemMixin + InfiniteSystemMixin diff --git a/kwant/_system.pyx b/kwant/_system.pyx index 71860ed6..2f58dd4f 100644 --- a/kwant/_system.pyx +++ b/kwant/_system.pyx @@ -1,4 +1,4 @@ -# Copyright 2011-2013 Kwant authors. +# Copyright 2011-2019 Kwant authors. # # This file is part of Kwant. It is subject to the license terms in the file # LICENSE.rst found in the top-level directory of this distribution and at @@ -12,12 +12,16 @@ import numpy as np from scipy import sparse as sp from itertools import chain import types +import bisect from .graph.core cimport CGraph, gintArraySlice from .graph.defs cimport gint from .graph.defs import gint_dtype from ._common import deprecate_args + +### Non-vectorized methods + msg = ('Hopping from site {0} to site {1} does not match the ' 'dimensions of onsite Hamiltonians of these sites.') @@ -352,3 +356,383 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None, mat = func(ham, args, params, self.graph, diag, from_sites, n_by_to_site, to_norb, to_off, from_norb, from_off) return (mat, to_norb, from_norb) if return_norb else mat + + +### Vectorized methods + + +_shape_error_msg = ( + "The following hoppings have matrix elements of incompatible shape " + "with other matrix elements: {}" +) + + +@cython.boundscheck(False) +def _vectorized_make_sparse(subgraphs, hams, long [:] norbs, long [:] orb_offsets, + long [:] site_offsets, shape=None): + ndata = sum(h.shape[0] * h.shape[1] * h.shape[2] for h in hams) + + cdef long [:] rows_view, cols_view + cdef complex [:] data_view + + rows = np.empty((ndata,), dtype=long) + cols = np.empty((ndata,), dtype=long) + data = np.empty((ndata,), dtype=complex) + rows_view = rows + cols_view = cols + data_view = data + + cdef long i, j, k, m, N, M, P, to_off, from_off,\ + ta, fa, to_norbs, from_norbs + cdef long [:] ts_offs, fs_offs + cdef complex [:, :, :] h + + m = 0 + # This outer loop zip() is pure Python, but that's ok, as it + # has very few entries and the inner loops are fully vectorized + for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams): + N = h.shape[0] + M = h.shape[1] + P = h.shape[2] + + if norbs[ta] != M or norbs[fa] != P: + to_sites = site_offsets[ta] + np.array(ts_offs) + from_sites = site_offsets[fa] + np.array(fs_offs) + hops = np.array([to_sites, from_sites]).transpose() + raise ValueError(_shape_error_msg.format(hops)) + + for i in range(N): + to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i] + from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i] + for j in range(M): + for k in range(P): + rows_view[m] = to_off + j + cols_view[m] = from_off + k + data_view[m] = h[i, j, k] + m += 1 + + if shape is None: + shape = (orb_offsets[-1], orb_offsets[-1]) + + return sp.coo_matrix((data, (rows, cols)), shape=shape) + + +@cython.boundscheck(False) +def _vectorized_make_dense(subgraphs, hams, long [:] norbs, long [:] orb_offsets, + long [:] site_offsets, shape=None): + if shape is None: + shape = (orb_offsets[-1], orb_offsets[-1]) + mat = np.zeros(shape, dtype=complex) + cdef complex [:, :] mat_view + mat_view = mat + + cdef long i, j, k, N, M, P, to_off, from_off,\ + ta, fa, to_norbs, from_norbs + cdef long [:] ts_offs, fs_offs + cdef complex [:, :, :] h + + # This outer loop zip() is pure Python, but that's ok, as it + # has very few entries and the inner loops are fully vectorized + for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams): + N = h.shape[0] + M = h.shape[1] + P = h.shape[2] + + if norbs[ta] != M or norbs[fa] != P: + to_sites = site_offsets[ta] + np.array(ts_offs) + from_sites = site_offsets[fa] + np.array(fs_offs) + hops = np.array([to_sites, from_sites]).transpose() + raise ValueError(_shape_error_msg.format(hops)) + + for i in range(N): + to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i] + from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i] + for j in range(M): + for k in range(P): + mat_view[to_off + j, from_off + k] = h[i, j, k] + return mat + + +def _count_norbs(syst, site_offsets, hams, args=(), params=None): + """Return the norbs and orbital offset of each site family in 'syst' + + Parameters + ---------- + syst : `kwant.system.System` + site_offsets : sequence of int + Contains the index of the first site for each site array + in 'syst.site_arrays'. + hams : sequence of arrays or 'None' + The Hamiltonian for each term in 'syst.terms'. 'None' + if the corresponding term has not been evaluated. + args, params : positional and keyword arguments to the system Hamiltonian + """ + site_ranges = syst.site_ranges + if site_ranges is None: + # NOTE: Remove this branch once we can rely on + # site families storing the norb information. + + site_arrays = syst.site_arrays + family_norbs = {s.family: None for s in site_arrays} + + # Compute the norbs per site family using already evaluated + # Hamiltonian terms where possible + for ham, t in zip(hams, syst.terms): + (to_w, from_w), _ = syst.subgraphs[t.subgraph] + if ham is not None: + family_norbs[site_arrays[to_w].family] = ham.shape[1] + family_norbs[site_arrays[from_w].family] = ham.shape[2] + + # Evaluate Hamiltonian terms where we do not already have them + for n, t in enumerate(syst.terms): + (to_w, from_w), _ = syst.subgraphs[t.subgraph] + to_fam = site_arrays[to_w].family + from_fam = site_arrays[from_w].family + if family_norbs[to_fam] is None or family_norbs[from_fam] is None: + ham = syst.hamiltonian_term(n, args=args, params=params) + family_norbs[to_fam] = ham.shape[1] + family_norbs[from_fam] = ham.shape[2] + + # This case is technically allowed by the format (some sites present + # but no hamiltonian terms that touch them) but is very unlikely. + if any(norbs is None for norbs in family_norbs.values()): + raise ValueError("Cannot determine the number of orbitals for " + "some site families.") + + orb_offset = 0 + site_ranges = [] + for start, site_array in zip(site_offsets, syst.site_arrays): + norbs = family_norbs[site_array.family] + site_ranges.append((start, norbs, orb_offset)) + orb_offset += len(site_array) * norbs + site_ranges.append((site_offsets[-1], 0, orb_offset)) + site_ranges = np.array(site_ranges) + + _, norbs, orb_offsets = site_ranges.transpose() + # The final (extra) element in orb_offsets is the total number of orbitals + assert len(orb_offsets) == len(syst.site_arrays) + 1 + + return norbs, orb_offsets + + +def _expand_norbs(compressed_norbs, site_offsets): + "Return norbs per site, given norbs per site array." + norbs = np.empty(site_offsets[-1], int) + for start, stop, norb in zip(site_offsets, site_offsets[1:], + compressed_norbs): + norbs[start:stop] = norb + return norbs + + +def _reverse_subgraph(subgraph): + (to_sa, from_sa), (to_off, from_off) = subgraph + return ((from_sa, to_sa), (from_off, to_off)) + + +@deprecate_args +@cython.binding(True) +def vectorized_hamiltonian_submatrix(self, args=(), sparse=False, + return_norb=False, *, params=None): + """Return The system Hamiltonian. + + Parameters + ---------- + args : tuple, defaults to empty + Positional arguments to pass to ``hamiltonian_term``. Mutually + exclusive with 'params'. + sparse : bool + Whether to return a sparse or a dense matrix. Defaults to ``False``. + return_norb : bool + Whether to return arrays of numbers of orbitals. Defaults to ``False``. + params : dict, optional + Dictionary of parameter names and their values. Mutually exclusive + with 'args'. + + Returns + ------- + hamiltonian_part : numpy.ndarray or scipy.sparse.coo_matrix + The Hamiltonian of the system. + norb : array of integers + Numbers of orbitals on each site. Only returned when ``return_norb`` + is true. + + Notes + ----- + Providing positional arguments via 'args' is deprecated, + instead, provide named parameters as a dictionary via 'params'. + """ + + site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays]) + + subgraphs = [ + self.subgraphs[t.subgraph] + for t in self.terms + ] + # Add Hermitian conjugates + subgraphs += [ + _reverse_subgraph(self.subgraphs[t.subgraph]) + for t in self.terms + if t.hermitian + ] + + hams = [ + self.hamiltonian_term(n, args=args, params=params) + for n, _ in enumerate(self.terms) + ] + # Add Hermitian conjugates + hams += [ + ham.conjugate().transpose(0, 2, 1) + for ham, t in zip(hams, self.terms) + if t.hermitian + ] + + norbs, orb_offsets = _count_norbs(self, site_offsets, hams, + args=args, params=params) + + func = _vectorized_make_sparse if sparse else _vectorized_make_dense + mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets) + + if return_norb: + return (mat, _expand_norbs(norbs, site_offsets)) + else: + return mat + + +@deprecate_args +@cython.binding(True) +def vectorized_cell_hamiltonian(self, args=(), sparse=False, *, params=None): + """Hamiltonian of a single cell of the infinite system. + + Providing positional arguments via 'args' is deprecated, + instead, provide named parameters as a dictionary via 'params'. + """ + + site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays]) + # Site array where next cell starts + next_cell = bisect.bisect(site_offsets, self.cell_size) - 1 + + def inside_cell(term): + (to_which, from_which), _= self.subgraphs[term.subgraph] + return to_which < next_cell and from_which < next_cell + + cell_terms = [ + n for n, t in enumerate(self.terms) + if inside_cell(t) + ] + + subgraphs = [ + self.subgraphs[self.terms[n].subgraph] + for n in cell_terms + ] + # Add Hermitian conjugates + subgraphs += [ + _reverse_subgraph(self.subgraphs[self.terms[n].subgraph]) + for n in cell_terms + if self.terms[n].hermitian + ] + + hams = [ + self.hamiltonian_term(n, args=args, params=params) + for n in cell_terms + ] + # Add Hermitian conjugates + hams += [ + ham.conjugate().transpose(0, 2, 1) + for ham, n in zip(hams, cell_terms) + if self.terms[n].hermitian + ] + + # _count_norbs requires passing hamiltonians for all terms, or 'None' + # if it has not been evaluated + all_hams = [None] * len(self.terms) + for n, ham in zip(cell_terms, hams): + all_hams[n] = ham + norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams, + args=args, params=params) + + shape = (orb_offsets[next_cell], orb_offsets[next_cell]) + func = _vectorized_make_sparse if sparse else _vectorized_make_dense + mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape) + + return mat + + +@deprecate_args +@cython.binding(True) +def vectorized_inter_cell_hopping(self, args=(), sparse=False, *, params=None): + """Hopping Hamiltonian between two cells of the infinite system. + + Providing positional arguments via 'args' is deprecated, + instead, provide named parameters as a dictionary via 'params'. + """ + + # Take advantage of the fact that there are distinct entries in + # onsite_terms for sites inside the unit cell, and distinct entries + # in onsite_terms for hoppings between the unit cell sites and + # interface sites. This way we only need to check the first entry + # in onsite/hopping_terms + + site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays]) + # Site array where next cell starts + next_cell = bisect.bisect(site_offsets, self.cell_size) - 1 + + inter_cell_hopping_terms = [ + n for n, t in enumerate(self.terms) + # *from* site is in interface + if (self.subgraphs[t.subgraph][0][1] >= next_cell + and self.subgraphs[t.subgraph][0][0] < next_cell) + ] + reversed_inter_cell_hopping_terms = [ + n for n, t in enumerate(self.terms) + # *to* site is in interface + if (self.subgraphs[t.subgraph][0][0] >= next_cell + and self.subgraphs[t.subgraph][0][1] < next_cell) + ] + + # Evaluate inter-cell hoppings only + inter_cell_hams = [ + self.hamiltonian_term(n, args=args, params=params) + for n in inter_cell_hopping_terms + ] + reversed_inter_cell_hams = [ + self.hamiltonian_term(n, args=args, params=params) + .conjugate().transpose(0, 2, 1) + for n in reversed_inter_cell_hopping_terms + ] + + hams = inter_cell_hams + reversed_inter_cell_hams + + subgraphs = [ + self.subgraphs[self.terms[n].subgraph] + for n in inter_cell_hopping_terms + ] + subgraphs += [ + _reverse_subgraph(self.subgraphs[self.terms[n].subgraph]) + for n in reversed_inter_cell_hopping_terms + ] + + # All the 'from' sites are in the previous domain, but to build a + # matrix we need to get the equivalent sites in the fundamental domain. + # We set the 'from' site array to the one from the fundamental domain. + subgraphs = [ + ((to_sa, from_sa - next_cell), (to_off, from_off)) + for (to_sa, from_sa), (to_off, from_off) in subgraphs + ] + + # _count_norbs requires passing hamiltonians for all terms, or 'None' + # if it has not been evaluated + all_hams = [None] * len(self.terms) + for n, ham in zip(inter_cell_hopping_terms, inter_cell_hams): + all_hams[n] = ham + for n, ham in zip(reversed_inter_cell_hopping_terms, + reversed_inter_cell_hams): + # Transpose to get back correct shape wrt. original term + all_hams[n] = ham.transpose(0, 2, 1) + norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams, + args=args, params=params) + + shape = (orb_offsets[next_cell], + orb_offsets[len(self.site_arrays) - next_cell]) + func = _vectorized_make_sparse if sparse else _vectorized_make_dense + mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape) + return mat diff --git a/kwant/system.py b/kwant/system.py index b08f27cf..a65f9259 100644 --- a/kwant/system.py +++ b/kwant/system.py @@ -10,7 +10,8 @@ __all__ = [ 'Site', 'SiteArray', 'SiteFamily', - 'System', 'FiniteSystem', 'InfiniteSystem', + 'System', 'VectorizedSystem', 'FiniteSystem', 'FiniteVectorizedSystem', + 'InfiniteSystem', 'InfiniteVectorizedSystem', ] import abc @@ -18,7 +19,7 @@ import warnings import operator from copy import copy from collections import namedtuple -from functools import total_ordering +from functools import total_ordering, lru_cache import numpy as np from . import _system from ._common import deprecate_args, KwantDeprecationWarning @@ -347,12 +348,100 @@ class System(metaclass=abc.ABCMeta): details = ', and '.join((', '.join(details[:-1]), details[-1])) return '<{} with {}>'.format(self.__class__.__name__, details) + hamiltonian_submatrix = _system.hamiltonian_submatrix -# Add a C-implemented function as an unbound method to class System. -System.hamiltonian_submatrix = _system.hamiltonian_submatrix +Term = namedtuple( + "Term", + ["subgraph", "hermitian", "parameters"], +) -class FiniteSystem(System, metaclass=abc.ABCMeta): + +class VectorizedSystem(System, metaclass=abc.ABCMeta): + """Abstract general low-level system with support for vectorization. + + Attributes + ---------- + graph : kwant.graph.CGraph + The system graph. + subgraphs : sequence of tuples + Each subgraph has the form '((idx1, idx2), (offsets1, offsets2))' + where 'offsets1' and 'offsets2' index sites within the site arrays + indexed by 'idx1' and 'idx2'. + terms : sequence of tuples + Each tuple has the following structure: + (subgraph: int, hermitian: bool, parameters: List(str)) + 'subgraph' indexes 'subgraphs' and supplies the to/from sites of this + term. 'hermitian' is 'True' if the term needs its Hermitian + conjugate to be added when evaluating the Hamiltonian, and 'parameters' + contains a list of parameter names used when evaluating this term. + site_arrays : sequence of SiteArray + The sites of the system. + site_ranges : None or Nx3 integer array + Has 1 row per site array, plus one extra row. Each row consists + of ``(first_site, norbs, orb_offset)``: the index of the first + site in the site array, the number of orbitals on each site in + the site array, and the offset of the first orbital of the first + site in the site array. In addition, the final row has the form + ``(len(graph.num_nodes), 0, tot_norbs)`` where ``tot_norbs`` is the + total number of orbitals in the system. ``None`` if any site array + in 'site_arrays' does not have 'norbs' specified. Note 'site_ranges' + is directly computable from 'site_arrays'. + parameters : frozenset of strings + The names of the parameters on which the system depends. This attribute + is provisional and may be changed in a future version of Kwant + + Notes + ----- + The sites of the system are indexed by integers ranging from 0 to + ``self.graph.num_nodes - 1``. + + Optionally, a class derived from ``System`` can provide a method ``pos`` which + is assumed to return the real-space position of a site given its index. + """ + @abc.abstractmethod + def hamiltonian_term(self, term_number, selector=slice(None), + args=(), params=None): + """Return the Hamiltonians for hamiltonian term number k. + + Parameters + ---------- + term_number : int + The number of the term to evaluate. + selector : slice or sequence of int, default: slice(None) + The elements of the term to evaluate. + args : tuple + Positional arguments to the term. (Deprecated) + params : dict + Keyword parameters to the term + + Returns + ------- + hamiltonian : 3d complex array + Has shape ``(N, P, Q)`` where ``N`` is the number of matrix + elements in this term (or the number selected by 'selector' + if provided), ``P`` and ``Q`` are the number of orbitals in the + 'to' and 'from' site arrays associated with this term. + + Providing positional arguments via 'args' is deprecated, + instead, provide named parameters as a dictionary via 'params'. + """ + @property + @lru_cache(1) + def site_ranges(self): + site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays]) + norbs = [arr.family.norbs for arr in self.site_arrays] + [0] + if any(norb is None for norb in norbs): + return None + orb_offsets = np.cumsum( + [0] + [len(arr) * arr.family.norbs for arr in self.site_arrays] + ) + return np.array([site_offsets, norbs, orb_offsets]).transpose() + + hamiltonian_submatrix = _system.vectorized_hamiltonian_submatrix + + +class FiniteSystemMixin(metaclass=abc.ABCMeta): """Abstract finite low-level system, possibly with leads. Attributes @@ -475,7 +564,15 @@ class FiniteSystem(System, metaclass=abc.ABCMeta): return symmetries.validate(ham) -class InfiniteSystem(System, metaclass=abc.ABCMeta): +class FiniteSystem(System, FiniteSystemMixin, metaclass=abc.ABCMeta): + pass + + +class FiniteVectorizedSystem(VectorizedSystem, FiniteSystemMixin, metaclass=abc.ABCMeta): + pass + + +class InfiniteSystemMixin(metaclass=abc.ABCMeta): """Abstract infinite low-level system. An infinite system consists of an infinite series of identical cells. @@ -516,30 +613,10 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta): infinite system. The other scheme has the numbers of site 0 and 1 exchanged, as well as of site 3 and 4. + Sites in the fundamental domain cell must belong to a different site array + than the sites in the previous cell. In the above example this means that + sites '(0, 1, 2)' and '(3, 4)' must belong to different site arrays. """ - @deprecate_args - def cell_hamiltonian(self, args=(), sparse=False, *, params=None): - """Hamiltonian of a single cell of the infinite system. - - Providing positional arguments via 'args' is deprecated, - instead, provide named parameters as a dictionary via 'params'. - """ - cell_sites = range(self.cell_size) - return self.hamiltonian_submatrix(args, cell_sites, cell_sites, - sparse=sparse, params=params) - - @deprecate_args - def inter_cell_hopping(self, args=(), sparse=False, *, params=None): - """Hopping Hamiltonian between two cells of the infinite system. - - Providing positional arguments via 'args' is deprecated, - instead, provide named parameters as a dictionary via 'params'. - """ - cell_sites = range(self.cell_size) - interface_sites = range(self.cell_size, self.graph.num_nodes) - return self.hamiltonian_submatrix(args, cell_sites, interface_sites, - sparse=sparse, params=params) - @deprecate_args def modes(self, energy=0, args=(), *, params=None): """Return mode decomposition of the lead @@ -623,6 +700,37 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta): return list(broken) +class InfiniteSystem(System, InfiniteSystemMixin, metaclass=abc.ABCMeta): + + @deprecate_args + def cell_hamiltonian(self, args=(), sparse=False, *, params=None): + """Hamiltonian of a single cell of the infinite system. + + Providing positional arguments via 'args' is deprecated, + instead, provide named parameters as a dictionary via 'params'. + """ + cell_sites = range(self.cell_size) + return self.hamiltonian_submatrix(args, cell_sites, cell_sites, + sparse=sparse, params=params) + + @deprecate_args + def inter_cell_hopping(self, args=(), sparse=False, *, params=None): + """Hopping Hamiltonian between two cells of the infinite system. + + Providing positional arguments via 'args' is deprecated, + instead, provide named parameters as a dictionary via 'params'. + """ + cell_sites = range(self.cell_size) + interface_sites = range(self.cell_size, self.graph.num_nodes) + return self.hamiltonian_submatrix(args, cell_sites, interface_sites, + sparse=sparse, params=params) + + +class InfiniteVectorizedSystem(VectorizedSystem, InfiniteSystemMixin, metaclass=abc.ABCMeta): + cell_hamiltonian = _system.vectorized_cell_hamiltonian + inter_cell_hopping = _system.vectorized_inter_cell_hopping + + class PrecalculatedLead: def __init__(self, modes=None, selfenergy=None): """A general lead defined by its self energy. -- GitLab