diff --git a/doc/source/pre/whatsnew/1.3.rst b/doc/source/pre/whatsnew/1.3.rst
index b5375c2f977fdfc5c15ca27cc5c3b3af6eba9465..996a9061bebf4a52c43f55238445983cdaf870bd 100644
--- a/doc/source/pre/whatsnew/1.3.rst
+++ b/doc/source/pre/whatsnew/1.3.rst
@@ -83,3 +83,9 @@ attach_lead() can now handle leads with greater than nearest-neighbor hoppings
When attaching a lead with greater than nearest-neighbor hoppings, the symmetry
period of the finalized lead is suitably extended and the unit cell size is
increased.
+
+Reference implementation of the Kernel Polynomial Method
+--------------------------------------------------------
+The kernel polynomial method is now implemented within Kwant to obtain the
+density of states or, more generally, the spectral density of a given operator
+acting on a system or Hamiltonian.
diff --git a/doc/source/reference/index.rst b/doc/source/reference/index.rst
index 6c2c2f9f3dec67c7dae310094aed044465500f9f..fd88a0a5db9b0f2322dbfc03cc26a38be19c5d41 100644
--- a/doc/source/reference/index.rst
+++ b/doc/source/reference/index.rst
@@ -37,3 +37,4 @@ The following modules provide functionality for special applications.
kwant.operator
kwant.digest
kwant.rmt
+ kwant.kpm
diff --git a/doc/source/reference/kwant.kpm.rst b/doc/source/reference/kwant.kpm.rst
new file mode 100644
index 0000000000000000000000000000000000000000..adec65d795005c7e1aaab2ac2931446f49cfeef6
--- /dev/null
+++ b/doc/source/reference/kwant.kpm.rst
@@ -0,0 +1,7 @@
+:mod:`kwant.kpm` -- Kernel Polynomial Method
+============================================
+
+.. automodule:: kwant.kpm
+ :members:
+ :special-members:
+ :exclude-members: __weakref__
diff --git a/kwant/__init__.py b/kwant/__init__.py
index a4ae67912c25193f721b8d0757d820eda0192e88..6dcc99eec2062e2befca40811efefa408e161c8e 100644
--- a/kwant/__init__.py
+++ b/kwant/__init__.py
@@ -30,7 +30,7 @@ __all__.extend(['KwantDeprecationWarning', 'UserCodeError'])
from ._common import version as __version__
for module in ['system', 'builder', 'lattice', 'solvers', 'digest', 'rmt',
- 'operator']:
+ 'operator', 'kpm']:
exec('from . import {0}'.format(module))
__all__.append(module)
diff --git a/kwant/kpm.py b/kwant/kpm.py
new file mode 100644
index 0000000000000000000000000000000000000000..3f7e661d8b7a9cc8d2013c28771539967cd4d041
--- /dev/null
+++ b/kwant/kpm.py
@@ -0,0 +1,534 @@
+# -*- coding: utf-8 -*-
+# Copyright 2011-2016 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
+# http://kwant-project.org/license. A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+__all__ = ['SpectralDensity']
+
+import warnings
+
+import numpy as np
+import scipy
+import scipy.sparse.linalg as sla
+import scipy.fftpack as fft
+
+from . import system
+from ._common import ensure_isinstance, ensure_rng
+
+
+class SpectralDensity:
+ """Calculate the spectral density of an operator.
+
+ The **SpectralDensity** class makes use of the kernel polynomial
+ method (KPM), presented in [1]_, to obtain the spectral density
+ :math:`Ï_A(e)`, as a function of the energy :math:`e`, of some
+ operator :math:`A` that acts on a kwant system or a Hamiltonian.
+ In general
+
+ .. math::
+ Ï_A(e) = Ï(e) A(e),
+
+ where :math:`Ï(e)` is the density of states, and :math:`A(e)` is
+ the expectation value of :math:`A` for all the eigenstates with
+ energy :math:`e`.
+
+ The **SpectralDensity** class can be called as
+
+ >>> rho = kwant.kpm.SpectralDensity(fsyst, operator=A),
+
+ with optional
+ arguments that specify the accuracy of the approximation. The
+ parameters ``num_moments``, ``num_rand_vecs``, and
+ ``num_sampling_points`` can be specified when creating an instance
+ or by calling
+
+ >>> rho.increase_accuracy(num_moments=200,
+ ... num_rand_vecs=10,
+ ... num_sampling_points=400)
+
+ Parameters
+ ----------
+ syst_or_ham : `~kwant.system.FiniteSystem` or matrix Hamiltonian
+ If a system is passed, it should contain no leads.
+ args : tuple, optional
+ Positional arguments to pass to the system.
+ operator : operator, dense matrix, or sparse matrix, optional
+ Operator for which the spectral density will be evaluated. If
+ it is callable, the ``densities`` at each energy will have the
+ dimension of the result of `operator(bra, ket)`. If it has a
+ ``dot`` method, such as ``numpy.ndarray`` and ``scipy.sparse.
+ matrices``, the densities will be scalars.
+ If no operator is provided, the density of states is calculated with a
+ faster algorithm.
+ num_rand_vecs : integer, default: 10
+ Number of random vectors for the KPM method.
+ num_moments : integer, default: 100
+ Number of moments, order of the KPM expansion.
+ num_sampling_points : integer, optional
+ Number of points where the spectral density will be evaluated.
+ If not provided, ``2\*num_moments`` will be used.
+ vector_factory : function, optional
+ The user defined function ``f(n)`` generates random vectors of
+ length ``n`` that will be used in the algorithm.
+ If not provided, random phase vectors are used.
+ The default random vectors are optimal for most cases, see the
+ discussions in [1]_ and [2]_.
+ bounds : tuple of floats (lower, upper), optional
+ Lower and upper bounds for all the eigenvalues of the system.
+ If provided, they save computation time.
+ rng : seed, or random number generator, optional
+ Random number generator used by ``vector_factory``.
+ If not provided, numpy's rng will be used; if it is an Integer,
+ it will be used to seed numpy's rng, and if it is a random
+ number generator, this is the one used.
+
+ Notes
+ -----
+ When passing an ``operator`` defined in `~kwant.operator`, the
+ result of ``operator(bra, ket)`` depends on the attribute ``sum``
+ of such operator. If ``sum=True``, densities will be scalars, that
+ is, total densities of the system. If ``sum=False`` the densities
+ will be arrays of the length of the system, that is, local
+ densities.
+
+ .. [1] `Rev. Mod. Phys., Vol. 78, No. 1 (2006)
+ <https://arxiv.org/abs/cond-mat/0504627>`_.
+ .. [2] `Phys. Rev. E 69, 057701 (2004)
+ <https://arxiv.org/abs/cond-mat/0401202>`_
+
+ Examples
+ --------
+ In the following example, we will obtain the density of states of a
+ graphene sheet, defined as a honeycomb lattice with first nearest
+ neighbors coupling.
+
+ We start by importing kwant and defining a
+ `~kwant.system.FiniteSystem`,
+
+ >>> import kwant
+ ...
+ >>> def circle(pos):
+ ... x, y = pos
+ ... return x**2 + y**2 < 100
+ ...
+ >>> lat = kwant.lattice.honeycomb()
+ >>> syst = kwant.Builder()
+ >>> syst[lat.shape(circle, (0, 0))] = 0
+ >>> syst[lat.neighbors()] = 1
+
+ and after finalizing the system, create an instance of
+ **SpectralDensity**
+
+ >>> fsyst = syst.finalized()
+ >>> rho = kwant.kpm.SpectralDensity(fsyst)
+
+ The ``energies`` and ``densities`` can be accesed with
+
+ >>> energies, densities = rho()
+
+ or
+
+ >>> energies, densities = rho.energies, rho.densities
+
+ Attributes
+ ----------
+ energies : array of floats
+ Array of sampling points with length ``num_sampling_points`` in
+ the range of the spectrum.
+ densities : array of floats
+ Spectral density of the ``operator`` evaluated at the energies.
+
+ """
+
+ def __init__(self, syst_or_ham, args=(), operator=None, num_rand_vecs=10,
+ num_moments=100, num_sampling_points=None,
+ vector_factory=None, bounds=None, rng=None):
+ rng = ensure_rng(rng)
+ # self.epsilon ensures that the rescaled Hamiltonian has a
+ # spectrum strictly in the interval (-1,1).
+ self.epsilon = None
+ # Check if syst_or_ham is a finalized System or a Hamiltonian
+ try:
+ self.ham = scipy.sparse.csr_matrix(syst_or_ham)
+ except:
+ try:
+ ensure_isinstance(syst_or_ham, system.System)
+ self.ham = scipy.sparse.csr_matrix(
+ syst_or_ham.hamiltonian_submatrix(args=args, sparse=True))
+ except TypeError:
+ raise ValueError('Parameter `syst_or_ham` is not a '
+ 'Hamiltonian neither a (finalized) '
+ '`kwant.system`.')
+ # Check if operator is a sparse matrix or None.
+ if operator is None:
+ self.operator = None
+ else:
+ if callable(operator):
+ self.operator = operator
+ elif hasattr(operator, 'dot') and not hasattr(operator, 'act'):
+ operator = scipy.sparse.csr_matrix(operator)
+ self.operator = (lambda bra, ket:
+ np.vdot(bra, operator.dot(ket)))
+ else:
+ raise ValueError('Parameter `operator` has no `.dot` '
+ 'attribute and is not callable.')
+ self.num_rand_vecs = num_rand_vecs
+ self.num_moments = num_moments
+ # Default number of sampling points
+ if num_sampling_points is None:
+ self.num_sampling_points = 2 * self.num_moments
+ elif num_sampling_points < self.num_moments:
+ raise ValueError('The number of sampling points should be larger '
+ 'than the number of moments.')
+ else:
+ self.num_sampling_points = num_sampling_points
+
+ # calculate raw moments
+ self._vector_factory = vector_factory or \
+ (lambda n: np.exp(2j * np.pi * rng.random_sample(n)))
+ # store this vector for reproducibility
+ self._v0 = self._vector_factory(self.ham.shape[0])
+ # spectrum bounds if provided
+ self._given_bounds = bounds
+ self._rand_vect_list = []
+ self._last_two_alphas = []
+ self._moments_list = [0.] * self.num_rand_vecs
+ self._update_moments_list(self.num_moments, self.num_rand_vecs)
+ # sum moments of all random vectors
+ moments = np.sum(np.asarray(self._moments_list), axis=0)
+ # divide by the number of random vectors
+ moments = moments / self.num_rand_vecs
+ # divide by the length of the vectors to normalize
+ moments = moments / self.ham.shape[0]
+ # divide by scale factor to reflect the integral rescaling
+ moments = moments / self._a
+
+ # obtain energies, densities, and gammas of Eq. 83
+ xk_rescaled, rho, self._gammas = _calc_fft_moments(
+ moments, self.num_sampling_points)
+ # energies to normal scale
+ self.energies = xk_rescaled*self._a + self._b
+ self.densities = rho
+
+ def __call__(self, energy=None):
+ """Return the spectral density evaluated at ``energy``.
+
+ Parameters
+ ----------
+ energy : float or sequence of float, optional
+
+ Returns
+ -------
+ float, if ``energy`` is float, array of float if ``energy``
+ is a sequence, a tuple (energies, densities) if
+ ``energy`` was not provided.
+
+ Notes
+ -----
+ If ``energy`` is not provided, then the densities are obtained by
+ Fast Fourier Transform of the Chebyshev moments.
+ """
+ if energy is None:
+ return self.energies, self.densities
+ else:
+ energy = np.asarray(energy, dtype=complex)
+ rescaled_energy = (energy - self._b) / self._a
+ g_e = np.pi * np.sqrt(1 - rescaled_energy) *\
+ np.sqrt(1 + rescaled_energy)
+
+ # calculate the coefficients for Chebyshev expansion
+ # sum moments of all random vectors
+ moments = np.sum(np.asarray(self._moments_list), axis=0)
+ # divide by the number of random vectors
+ moments = moments / self.num_rand_vecs
+ # divide by the length of the vectors to normalize
+ moments = moments / self.ham.shape[0]
+ # divide by scale factor to reflect the integral rescaling
+ moments = moments / self._a
+
+ m = np.arange(self.num_moments)
+
+ kernel = ((self.num_moments - m + 1) *
+ np.cos(np.pi * m/(self.num_moments + 1)) +
+ np.sin(np.pi * m/(self.num_moments + 1)) /
+ np.tan(np.pi/(self.num_moments + 1)))
+ kernel = kernel / (self.num_moments + 1)
+
+ coef_cheb = np.zeros_like(moments)
+ coef_cheb[0] = moments[0]
+ coef_cheb[1:] = 2 * moments[1:] * kernel[1:]
+
+ return (np.polynomial.chebyshev.chebval(
+ rescaled_energy, coef_cheb) / g_e).real
+
+ def average(self, distribution_function=None):
+ """Returns the total spectral density.
+
+ Returns the integral over the whole spectrum with an optional
+ distribution function. ``distribution_function`` should be able to
+ take arrays as input. Defined using Gauss-Chebyshev
+ integration.
+ """
+ # This factor divides the sum to normalize the Gauss integral
+ # and rescales the integral back with ``self._a`` to normal scale.
+ factor = self._a/self.num_sampling_points
+ if distribution_function is None:
+ rho = self._gammas
+ else:
+ # The evaluation of the distribution function should be at
+ # the energies without rescaling.
+ distribution_array = distribution_function(self.energies)
+ rho = np.transpose(self._gammas.transpose() * distribution_array)
+ return factor * np.sum(rho, axis=0)
+
+ def increase_energy_resolution(self, tol=0.05, increase_num_moments=True):
+ """Set the minimum energy resolution to ``tol``.
+
+ The energy resolution is increased by increasing the number of
+ sampling points. By default the number of moments is also increased.
+ """
+ resolution = 2 * self._a / (self.num_sampling_points / 1.6)
+ if tol > resolution:
+ warnings.warn('Energy resolution is already smaller than tol.')
+ else:
+ num_sampling_points = np.ceil((1.6 * 2*self._a) / tol)
+ num_sampling_points = num_sampling_points.astype(int)
+ if increase_num_moments:
+ num_moments = num_sampling_points // 2
+ self.increase_accuracy(num_moments=num_moments,
+ num_sampling_points=num_sampling_points)
+ else:
+ self.increase_accuracy(num_sampling_points=num_sampling_points)
+
+ def increase_accuracy(self, num_moments=None, num_rand_vecs=None,
+ num_sampling_points=None):
+ """Increase the number of moments, random vectors, or sampling
+ points.
+
+ Parameters
+ ----------
+ num_moments, num_rand_vecs, num_sampling_points : integer
+ Number of Chebyshev moments, random vectors used for
+ sampling, and sampling points for the calculation of
+ densities.
+ """
+ # Calls to update_moments_list should be in this order:
+ # 1st. update random vectors and recalculate current moments,
+ # 2nd. update moments with the updated random vectors,
+ # 3rd. update sampling points and FFT densities.
+ num_rand_vecs = num_rand_vecs or self.num_rand_vecs
+ self._update_moments_list(self.num_moments, num_rand_vecs)
+ self.num_rand_vecs = num_rand_vecs
+
+ num_moments = num_moments or self.num_moments
+ self._update_moments_list(num_moments, self.num_rand_vecs)
+ self.num_moments = num_moments
+
+ num_sampling_points = num_sampling_points or self.num_sampling_points
+ if num_sampling_points < self.num_moments:
+ raise ValueError(
+ 'The number of sampling points should be larger than the '
+ 'number of moments.')
+ self.num_sampling_points = num_sampling_points
+
+ # sum moments of all random vectors
+ moments = np.sum(np.asarray(self._moments_list), axis=0)
+ # divide by the number of random vectors
+ moments = moments / self.num_rand_vecs
+ # divide by the length of the vectors to normalize
+ moments = moments / self.ham.shape[0]
+ # divide by scale factor to reflect the integral rescaling
+ moments = moments / self._a
+
+ xk_rescaled, rho, self._gammas = _calc_fft_moments(
+ moments, self.num_sampling_points)
+
+ # Sampling points to normal scale
+ self.energies = xk_rescaled * self._a + self._b
+ self.densities = rho
+
+ def _update_moments_list(self, n_moments, n_rand):
+ """Calculate the Chebyshev moments of an operator's spectral
+ density.
+
+ The algorithm is based on the KPM method as depicted in `Rev.
+ Mod. Phys., Vol. 78, No. 1 (2006)
+ <https://arxiv.org/abs/cond-mat/0504627>`_.
+
+ Parameters
+ ----------
+ n_moments : integer
+ Number of Chebyshev moments.
+ n_rand : integer
+ Number of random vectors used for sampling.
+ """
+ # Hamiltonian rescaled as in Eq. (24)
+ H_rescaled, self._a, self._b = _rescale(self.ham,
+ epsilon=self.epsilon,
+ v0=self._v0,
+ bounds=self._given_bounds)
+ self.bounds = (self._b - self._a, self._b + self._a)
+ dim = H_rescaled.shape[0]
+
+ if len(self._rand_vect_list) == 0:
+ for r in range(n_rand):
+ self._rand_vect_list.append(self._vector_factory(dim))
+ r_start = 0
+ r_stop = self.num_rand_vecs
+ else:
+ if self.num_rand_vecs == n_rand:
+ r_start = 0
+ r_stop = self.num_rand_vecs
+ elif self.num_rand_vecs < n_rand:
+ new_rand_vect = n_rand - self.num_rand_vecs
+ for r in range(new_rand_vect):
+ self._rand_vect_list.append(self._vector_factory(dim))
+ self._moments_list.extend([0.] * new_rand_vect)
+ r_start = self.num_rand_vecs
+ r_stop = n_rand
+ else:
+ warnings.warn('Decreasing number of random vectors')
+ r_start = 0
+ r_stop = n_rand
+
+ if n_moments == self.num_moments:
+ start = 2
+ new_moments = 0
+ else:
+ start = self.num_moments
+ new_moments = n_moments - self.num_moments
+ stop = n_moments
+
+ for r in range(r_start, r_stop):
+ alpha_zero = self._rand_vect_list[r]
+
+ one_moment = [0.] * n_moments
+ if new_moments == 0:
+ alpha = alpha_zero
+ alpha_next = H_rescaled.matvec(alpha)
+ if self.operator is None:
+ one_moment[0] = np.vdot(alpha_zero, alpha_zero)
+ one_moment[1] = np.vdot(alpha_zero, alpha_next)
+ else:
+ one_moment[0] = self.operator(alpha_zero, alpha_zero)
+ one_moment[1] = self.operator(alpha_zero, alpha_next)
+
+ elif new_moments > 0:
+ (alpha, alpha_next) = self._last_two_alphas[r]
+ one_moment[0:self.num_moments] = self._moments_list[r]
+ else:
+ raise ValueError('Attempt to decrease the number of moments.')
+ # Iteration over the moments
+ # Two cases can occur, depicted in Eq. (28) and in Eq. (29),
+ # respectively.
+ # ----
+ # In the first case, self.operator is None and we can use
+ # Eqs. (34) and (35) to obtain the density of states, with
+ # two moments ``one_moment`` for every new alpha.
+ # ----
+ # In the second case, the operator is not None and a matrix
+ # multiplication should be used.
+ if self.operator is None:
+ for n in range(start//2, stop//2):
+ alpha_save = alpha_next
+ alpha_next = 2 * H_rescaled.matvec(alpha_next) - alpha
+ alpha = alpha_save
+ # Following Eqs. (34) and (35)
+ one_moment[2*n] = 2 * np.vdot(alpha, alpha) -\
+ one_moment[0]
+ one_moment[2*n+1] = 2 * np.vdot(alpha_next, alpha) -\
+ one_moment[1]
+ else:
+ for n in range(start, stop):
+ alpha_save = alpha_next
+ alpha_next = 2 * H_rescaled.matvec(alpha_next) - alpha
+ alpha = alpha_save
+ one_moment[n] = self.operator(alpha_zero, alpha_next)
+
+ self._last_two_alphas.append((alpha, alpha_next))
+ self._moments_list[r] = np.asarray(one_moment).real
+
+
+# ### Auxiliary functions
+
+
+def _rescale(ham, epsilon=None, v0=None, bounds=None):
+ """Rescale a Hamiltonian and return a LinearOperator
+
+ Parameters
+ ----------
+ ham : 2D array
+ Hamiltonian of the system.
+ epsilon : scalar
+ Ensures that the bounds 'a' and 'b' are strict.
+ v0 : random vector
+ Used as the initial residual vector for the algorithm that
+ finds the lowest and highest eigenvalues.
+ bounds : tuple or None
+ Boundaries of the spectrum. If not provided the maximum and
+ minimum eigenvalues are calculated.
+ """
+ # calculate upper and lower bounds of the spectrum
+ TOL = 1e-14 # numerical precision
+ TOL_SP = 1e-6
+
+ epsilon = epsilon or 0.05
+ if bounds:
+ lmin = bounds[0]
+ lmax = bounds[1]
+ else:
+ lmax = float(sla.eigsh(ham, k=1, which='LA', return_eigenvectors=False,
+ tol=TOL_SP, v0=v0))
+ lmin = float(sla.eigsh(ham, k=1, which='SA', return_eigenvectors=False,
+ tol=TOL_SP, v0=v0))
+ a = np.abs(lmax-lmin) / (2. - epsilon)
+ b = (lmax+lmin) / 2.
+
+ if a < TOL:
+ raise ValueError(
+ 'The Hamiltonian has a single eigenvalue, it is not possible to '
+ 'obtain a spectral density.')
+
+ def rescaled(v):
+ return (ham.dot(v) - b * v) / a
+
+ return sla.LinearOperator(shape=ham.shape, matvec=rescaled), a, b
+
+
+def _calc_fft_moments(moments, n_sampling):
+ """This function takes the normalised moments and returns an array
+ of points and an array of the evaluated function at those points.
+ """
+ moments = np.asarray(moments)
+ # extra_shape handles the case where operators have vector outputs
+ n_moments, *extra_shape = moments.shape
+ moments_ext = np.zeros([n_sampling] + extra_shape)
+
+ # Jackson kernel, as in Eq. (71), and kernel improved moments,
+ # as in Eq. (81).
+ m = np.arange(n_moments)
+ kernel = ((n_moments - m + 1) * np.cos(np.pi * m / (n_moments + 1)) +
+ np.sin(np.pi * m / (n_moments + 1)) /
+ np.tan(np.pi / (n_moments + 1))) / (n_moments + 1)
+
+ # special points at the abscissas of Chebyshev integration
+ k = np.arange(0, n_sampling)
+ xk_rescaled = np.cos(np.pi * (k + 0.5) / n_sampling)
+ # prefactor in Eq. (80)
+ gk = np.pi * np.sqrt(1 - xk_rescaled ** 2)
+
+ # transposes handle the case where operators have vector outputs
+ moments_ext[:n_moments] = np.transpose(moments.transpose() * kernel)
+ # The function evaluated in these special data points is the FFT of
+ # the moments as in Eq. (83).
+ gammas = np.transpose(fft.dct(moments_ext.transpose(), type=3))
+
+ # Element-wise division of moments_FFT over gk, as in Eq. (83).
+ rho = np.transpose(np.divide(gammas.transpose(), gk))
+
+ # Reverse energies and densities to set ascending order.
+ return xk_rescaled[::-1], rho[::-1], gammas[::-1]
diff --git a/kwant/tests/test_kpm.py b/kwant/tests/test_kpm.py
new file mode 100644
index 0000000000000000000000000000000000000000..63c42471975316011c50ba9a164c305d4017f1dc
--- /dev/null
+++ b/kwant/tests/test_kpm.py
@@ -0,0 +1,522 @@
+# Copyright 2011-2016 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
+# http://kwant-project.org/license. A list of Kwant authors can be found in
+# the file AUTHORS.rst at the top-level directory of this distribution and at
+# http://kwant-project.org/authors.
+
+from copy import copy as copy
+from types import SimpleNamespace
+
+import pytest
+import numpy as np
+import scipy.sparse.linalg as sla
+from scipy.integrate import simps
+
+import kwant
+from ..kpm import _rescale
+from .._common import ensure_rng
+
+SpectralDensity = kwant.kpm.SpectralDensity
+
+# ### General parameters
+dim = 20
+p = SimpleNamespace(
+ num_moments=200,
+ num_rand_vecs=5,
+ num_sampling_points=400
+ )
+ham = kwant.rmt.gaussian(dim)
+
+# ### Auxiliary functions
+TOL = 1e-12
+TOL_SP = 1e-6
+TOL_WEAK = 1e-3
+
+
+def assert_allclose(arr1, arr2):
+ np.testing.assert_allclose(arr1, arr2, rtol=0., atol=TOL)
+
+
+def assert_allclose_sp(arr1, arr2):
+ np.testing.assert_allclose(arr1, arr2, rtol=0., atol=TOL_SP)
+
+
+def make_spectrum(ham, p, operator=None, vector_factory=None, rng=None):
+ """Create an instance of SpectralDensity class."""
+ return SpectralDensity(
+ ham,
+ operator=operator,
+ vector_factory=vector_factory,
+ num_moments=p.num_moments,
+ num_rand_vecs=p.num_rand_vecs,
+ num_sampling_points=p.num_sampling_points,
+ rng=rng
+ )
+
+
+def make_chain(r=dim, t=-1):
+ syst = kwant.Builder()
+ lat = kwant.lattice.chain(norbs=1)
+ for i in range(r):
+ syst[lat(i)] = 0
+ syst[lat.neighbors()] = t
+ return syst.finalized()
+
+
+def kpm_derivative(spectrum, e, order=1):
+ """Calculate the coeficients for Chebyshev expansion and its derivates."""
+ e = np.asarray(e, dtype=complex)
+ rescaled_energy = (e - spectrum._b) / spectrum._a
+ g_e = np.pi * np.sqrt(1 - rescaled_energy) * np.sqrt(1 + rescaled_energy)
+
+ m = np.arange(spectrum.num_moments)
+ kernel = ((spectrum.num_moments - m + 1) *
+ np.cos(np.pi * m / (spectrum.num_moments + 1)) +
+ np.sin(np.pi * m / (spectrum.num_moments + 1)) /
+ np.tan(np.pi / (spectrum.num_moments + 1)))
+ kernel = kernel / (spectrum.num_moments + 1)
+
+ moments = np.sum(spectrum._moments_list, axis=0) /\
+ (spectrum.num_rand_vecs * spectrum.ham.shape[0])
+ coef_cheb = np.zeros_like(moments)
+ coef_cheb[0] = moments[0]
+ coef_cheb[1:] = 2 * moments[1:] * kernel[1:]
+
+ i = 1
+ while i <= order:
+ i += 1
+ coef_cheb = np.polynomial.chebyshev.chebder(coef_cheb)
+
+ return np.real(np.polynomial.chebyshev.chebval(
+ rescaled_energy, coef_cheb)/g_e)
+
+
+def make_spectrum_and_peaks(ham, precise, threshold=0.005):
+ """Calculate the spectrum and the peaks in the array energies."""
+ spectrum = make_spectrum(ham, precise)
+ derivate_array = kpm_derivative(spectrum, spectrum.energies)
+ second_derivate_array = kpm_derivative(spectrum,
+ spectrum.energies, order=2)
+ # where the derivative changes sign
+ zero_index = np.where(np.diff(np.sign(derivate_array)))[0]
+ scale = threshold * np.max(spectrum.densities)
+ # where the spectrum is avobe the threshold
+ filter_index = zero_index[
+ np.where(spectrum.densities[zero_index] > scale)[0]]
+ # where is a local maxmimum
+ filter_index2 = filter_index[
+ np.where(second_derivate_array[filter_index] < 0)[0]]
+ return spectrum, filter_index2
+
+
+def deviation_from_eigenvalues(dim, precise):
+ """Return the maxmimum difference between the peaks in the spectrum and
+ the eigenvalues."""
+ ham = kwant.rmt.gaussian(dim)
+ spectrum, filter_index = make_spectrum_and_peaks(ham, precise)
+ try:
+ return 1./(2 * spectrum._a) * np.max(np.abs(
+ spectrum.energies[filter_index] - np.linalg.eigh(ham)[0]))
+ except ValueError:
+ """This means that two eigenvalues are nearly degenerate, therefore the
+ maximum deviation is the minimum difference between a peak and all the
+ eigenvalues."""
+ return np.max([1./(2*spectrum._a) *
+ np.min(np.abs(spectrum.energies[filter_index] -
+ np.linalg.eigh(ham)[0][i]))
+ for i in range(len(ham))])
+
+
+def find_peaks_in_an_array(array, threshold=0.005):
+ """Find the peaks above threshold and return the index inside array."""
+ derivate_array = np.gradient(array)
+ second_derivate_array = np.gradient(derivate_array)
+ # where the derivative changes sign
+ zero_index = np.where(np.diff(np.sign(derivate_array)))[0]
+ scale = threshold * np.max(array)
+ # where the spectrum is avobe the threshold
+ filter_index = zero_index[np.where(array[zero_index] > scale)[0]]
+ # where is a local maxmimum
+ filter_index2 = filter_index[
+ np.where(second_derivate_array[filter_index] < 0)[0]]
+ return filter_index2
+
+
+# ### Tests for consistency in top-level API ###
+
+
+def test_api_ham():
+ not_a_hamiltonian = str('not a hamiltonian')
+ with pytest.raises(ValueError):
+ SpectralDensity(not_a_hamiltonian)
+
+
+def test_api_operator():
+ not_an_operator = str('not an operator')
+ with pytest.raises(ValueError):
+ SpectralDensity(kwant.rmt.gaussian(dim),
+ operator=not_an_operator)
+
+
+def test_api_lower_num_moments():
+ spectrum = SpectralDensity(kwant.rmt.gaussian(dim))
+ num_moments = spectrum.num_moments - 1
+ with pytest.raises(ValueError):
+ spectrum.increase_accuracy(num_moments=num_moments)
+
+
+def test_api_lower_sampling_points():
+ spectrum = SpectralDensity(kwant.rmt.gaussian(dim))
+ num_samplint_points = spectrum.num_moments - 1
+ with pytest.raises(ValueError):
+ spectrum.increase_accuracy(num_sampling_points=num_samplint_points)
+
+
+def test_api_warning_less_sampling_points():
+ precise = copy(p)
+ precise.num_sampling_points = precise.num_moments - 1
+ ham = kwant.rmt.gaussian(dim)
+ with pytest.raises(ValueError):
+ make_spectrum(ham, precise)
+
+
+def test_api_lower_energy_resolution():
+ spectrum = make_spectrum(ham, p)
+ tol = np.max(np.abs(np.diff(spectrum.energies))) * 2.
+ with pytest.warns(UserWarning):
+ spectrum.increase_energy_resolution(tol=tol)
+
+
+def test_api_lower_num_rand_vecs():
+ spectrum = SpectralDensity(kwant.rmt.gaussian(dim))
+ num_rand_vecs = spectrum.num_rand_vecs - 1
+ with pytest.warns(UserWarning):
+ spectrum.increase_accuracy(num_rand_vecs=num_rand_vecs)
+
+
+def test_api_single_eigenvalue_error():
+ with pytest.raises(ValueError):
+ SpectralDensity(np.identity(dim, dtype=complex))
+
+
+def test_operator_none():
+ """Check operator=None gives the same results as operator=np.identity(),
+ with the same random vectors.
+ """
+ ham = kwant.rmt.gaussian(dim)
+ identity = np.identity(dim)
+
+ sp1 = SpectralDensity(ham, operator=None, rng=1)
+ sp2 = SpectralDensity(ham, operator=identity, rng=1)
+
+ # different algorithms are used so these arrays are equal up to TOL_SP
+ assert_allclose_sp(sp1.densities, sp2.densities)
+
+
+def test_bounds():
+ """Check operator=None gives the same results as operator=np.identity(),
+ with the same random vectors.
+ """
+ ham = kwant.rmt.gaussian(dim)
+ TOL_SP = 1e-6
+ rng = ensure_rng(1)
+ sp1 = SpectralDensity(ham, rng=rng)
+ # re initialize to obtain the same vector v0
+ rng = ensure_rng(1)
+ v0 = np.exp(2j * np.pi * rng.random_sample(dim))
+ lmax = float(sla.eigsh(
+ ham, k=1, which='LA', return_eigenvectors=False, tol=TOL_SP, v0=v0))
+ lmin = float(sla.eigsh(
+ ham, k=1, which='SA', return_eigenvectors=False, tol=TOL_SP, v0=v0))
+ sp2 = SpectralDensity(ham, bounds=(lmin, lmax), rng=1)
+
+ # different algorithms are used so these arrays are equal up to TOL_SP
+ assert_allclose_sp(sp1.densities, sp2.densities)
+
+def test_operator_user():
+ """Check operator=None gives the same results as operator=np.identity(),
+ with the same random vectors.
+ """
+ ham = kwant.rmt.gaussian(dim)
+ identity = np.identity(dim)
+ user_op = lambda bra, ket: np.vdot(bra, np.dot(identity, ket))
+
+ sp1 = SpectralDensity(ham, operator=user_op, rng=1)
+ sp2 = SpectralDensity(ham, operator=identity, rng=1)
+
+ # different algorithms are used so these arrays are equal up to TOL_SP
+ assert_allclose_sp(sp1.densities, sp2.densities)
+
+
+def test_kwant_syst():
+ """Check that when a kwant system is passed, the results are the same as
+ for the Hamiltonian (in both sparse and dense formats) of the same kwant
+ system.
+ """
+ syst = make_chain()
+ spectrum_syst = make_spectrum(syst, p, rng=1)
+ ham_sparse = syst.hamiltonian_submatrix(sparse=True)
+ ham_dense = syst.hamiltonian_submatrix(sparse=False)
+
+ spectrum_sparse = make_spectrum(ham_sparse, p, rng=1)
+ spectrum_dense = make_spectrum(ham_dense, p, rng=1)
+
+ # same algorithms are used so these arrays are equal up to TOL
+ assert_allclose(spectrum_syst.densities, spectrum_sparse.densities)
+ assert_allclose(spectrum_syst.densities, spectrum_dense.densities)
+
+
+def test_kwant_op():
+ """Check that the kwant.operator.Density gives the same result as the
+ identity operator.
+ """
+ syst = make_chain()
+ kwant_op = kwant.operator.Density(syst)
+ spectrum_syst = make_spectrum(syst, p, operator=kwant_op, rng=1)
+
+ ham = syst.hamiltonian_submatrix()
+ identity = np.identity(dim)
+ spectrum = make_spectrum(ham, p, operator=identity, rng=1)
+
+ assert spectrum_syst.densities.shape[1] == ham.shape[0]
+ # same algorithms are used so these arrays are equal up to TOL
+ assert_allclose(np.sum(spectrum_syst.densities, axis=1),
+ spectrum.densities)
+
+def test_kwant_op_average():
+ """Check that the kwant.operator.Density gives the same result as the
+ identity operator.
+ """
+ syst = make_chain()
+ kwant_op = kwant.operator.Density(syst)
+ spectrum_syst = make_spectrum(syst, p, operator=kwant_op, rng=1)
+
+ ham = syst.hamiltonian_submatrix()
+ identity = np.identity(dim)
+ spectrum = make_spectrum(ham, p, operator=identity, rng=1)
+ ones = lambda x: np.ones_like(x)
+
+ assert spectrum_syst.densities.shape[1] == ham.shape[0]
+ # same algorithms are used so these arrays are equal up to TOL
+ assert_allclose(np.sum(spectrum_syst.average(distribution_function=ones)),
+ spectrum.average())
+
+
+# ## test for methods to work as expected
+
+# ### increase num_moments
+
+
+def test_increase_num_moments():
+ spectrum = make_spectrum(ham, p, rng=1)
+ precise = copy(p)
+ precise.num_moments = 2 * p.num_moments
+
+ spectrum_raise = make_spectrum(ham, precise, rng=1)
+ spectrum.increase_accuracy(num_moments=precise.num_moments)
+
+ # Check bit for bit equality
+ assert np.all(spectrum_raise.densities == spectrum.densities)
+
+# ### increase num_moments with an operator (different algorithm)
+
+
+def test_increase_num_moments_op():
+ ham = kwant.rmt.gaussian(dim)
+ identity = np.identity(dim)
+
+ sp1 = make_spectrum(ham, p, operator=None, rng=1)
+ sp2 = make_spectrum(ham, p, operator=identity, rng=1)
+
+ sp1.increase_accuracy(num_moments=2*sp1.num_moments)
+ sp2.increase_accuracy(num_moments=2*sp2.num_moments)
+
+ # different algorithms are used so these arrays are equal up to TOL_SP
+ assert_allclose_sp(sp1.densities, sp2.densities)
+
+# ### increase num_random_vecs
+
+
+def test_increase_num_rand_vecs():
+ precise = copy(p)
+ precise.num_rand_vecs = 2 * p.num_rand_vecs
+
+ spectrum_raise = make_spectrum(ham, precise, rng=1)
+ spectrum = make_spectrum(ham, p, rng=1)
+ spectrum.increase_accuracy(num_rand_vecs=precise.num_rand_vecs)
+ # Check bit for bit equality
+
+ assert np.all(spectrum_raise.densities == spectrum.densities)
+
+# ### increase num_sampling_points
+
+
+def test_increase_num_sampling_points():
+ precise = copy(p)
+ precise.sampling_points = 2 * p.num_sampling_points
+
+ spectrum_raise = make_spectrum(ham, precise, rng=1)
+ spectrum = make_spectrum(ham, p, rng=1)
+ spectrum.increase_accuracy(num_moments=precise.num_moments)
+
+ # Check bit for bit equality
+ assert np.all(spectrum_raise.densities == spectrum.densities)
+
+# ### Convergence for higher number of moments ###
+
+
+def test_check_convergence_decreasing_values():
+ difference_list = []
+ depth = 2
+
+ for i in range(depth):
+ precise = SimpleNamespace(
+ num_moments=10*dim + 100*i*dim,
+ num_rand_vecs=dim//2,
+ num_sampling_points=None
+ )
+ results = []
+ np.random.seed(1)
+ iterations = 3
+
+ for ii in range(iterations):
+ results.append(deviation_from_eigenvalues(dim, precise))
+
+ difference_list.append(results)
+
+ assert(np.all(np.diff(np.sum(difference_list, axis=1) / iterations) < 0.))
+ assert(np.sum(difference_list, axis=1)[-1] / iterations < TOL_WEAK)
+
+# ### Convergence for a custom vector_factory
+
+
+def test_convergence_custom_vector_factory():
+ rng = ensure_rng(1)
+
+ def random_binary_vectors(dim):
+ return np.rint(rng.random_sample(dim)) * 2 - 1
+
+ # extracted from `deviation_from_eigenvalues
+ def deviation(ham, spectrum):
+ filter_index = find_peaks_in_an_array(spectrum.densities)
+ try:
+ return 1./(2 * spectrum._a) * np.max(np.abs(
+ spectrum.energies[filter_index] - np.linalg.eigh(ham)[0]))
+ except ValueError:
+ """This means that two eigenvalues are nearly degenerate, therefore
+ the maximum deviation is the minimum difference between a peak and
+ all the eigenvalues."""
+
+ return np.max([1./(2 * spectrum._a) *
+ np.min(np.abs(spectrum.energies[filter_index] -
+ np.linalg.eigh(ham)[0][i]))
+ for i in range(len(ham))])
+
+ difference_list = []
+ depth = 2
+
+ for i in range(depth):
+ precise = SimpleNamespace(
+ num_moments=10*dim + 100*i*dim,
+ num_rand_vecs=dim//2,
+ num_sampling_points=None
+ )
+ results = []
+ iterations = 3
+
+ for ii in range(iterations):
+ ham = kwant.rmt.gaussian(dim)
+ spectrum = make_spectrum(ham, precise,
+ vector_factory=random_binary_vectors)
+ results.append(deviation(ham, spectrum))
+
+ difference_list.append(results)
+
+ assert np.all(np.diff(np.sum(difference_list, axis=1) / iterations) < 0.)
+ assert np.sum(difference_list, axis=1)[-1]/iterations < TOL_WEAK
+
+# ## Consistency checks for internal functions
+
+# ### check call
+
+
+def test_call_no_argument():
+ ham = kwant.rmt.gaussian(dim)
+ spectrum = make_spectrum(ham, p)
+ energies, densities = spectrum()
+
+ # different algorithms are used so these arrays are equal up to TOL_SP
+ assert_allclose_sp(energies, spectrum.energies)
+ assert_allclose_sp(densities, spectrum.densities)
+
+
+def test_call():
+ ham = kwant.rmt.gaussian(dim)
+ spectrum = make_spectrum(ham, p)
+ densities_array = np.array([spectrum(e) for e in spectrum.energies])
+
+ # different algorithms are used so these arrays are equal up to TOL_SP
+ assert_allclose_sp(densities_array, spectrum.densities)
+
+# ### check average
+
+
+def test_average():
+ ham = kwant.rmt.gaussian(dim)
+ spectrum = make_spectrum(ham, p)
+ ones = lambda x: np.ones_like(x)
+ assert np.abs(spectrum.average() - simps(spectrum.densities,
+ x=spectrum.energies)) < TOL_SP
+ assert np.abs(spectrum.average() - spectrum.average(
+ distribution_function=ones)) < TOL
+
+# ### check increase_energy_resolution
+
+
+def test_increase_energy_resolution():
+ spectrum, filter_index = make_spectrum_and_peaks(ham, p)
+
+ old_sampling_points = spectrum.num_sampling_points
+ tol = np.max(np.abs(np.diff(spectrum.energies))) / 2
+
+ spectrum.increase_energy_resolution(tol=tol)
+ new_sampling_points = spectrum.num_sampling_points
+
+ assert old_sampling_points < new_sampling_points
+ assert np.max(np.abs(np.diff(spectrum.energies))) < tol
+
+
+def test_increase_energy_resolution_no_moments():
+ spectrum, filter_index = make_spectrum_and_peaks(ham, p)
+
+ old_sampling_points = spectrum.num_sampling_points
+ tol = np.max(np.abs(np.diff(spectrum.energies))) / 2
+
+ spectrum.increase_energy_resolution(tol=tol, increase_num_moments=False)
+ new_sampling_points = spectrum.num_sampling_points
+
+ assert old_sampling_points < new_sampling_points
+ assert np.max(np.abs(np.diff(spectrum.energies))) < tol
+
+# ### check _rescale
+
+
+def test_rescale():
+ ham = kwant.rmt.gaussian(dim)
+ spectrum, filter_index = make_spectrum_and_peaks(ham, p)
+
+ ham_operator, a, b = _rescale(ham)
+ rescaled_eigvalues, rescaled_eigvectors = np.linalg.eigh((
+ ham - b * np.identity(len(ham))) / a)
+
+ assert np.all(-1 < rescaled_eigvalues) and np.all(rescaled_eigvalues < 1)
+
+ """The test
+ `assert np.max(np.abs(eigvectors - rescaled_eigvectors)) < TOL`
+ is not good because eigenvectors could be defined up to a phase. It should
+ compare that the eigenvectors are proportional to eachother. One way to
+ test this is that the product gives a complex number in the unit circle."""
+ eigvalues, eigvectors = np.linalg.eigh(ham)
+ assert np.all(1 - np.abs(np.vdot(eigvectors, rescaled_eigvectors)) < TOL)