From 3f675e337cd28a5216d7c1c682aa19b6ffe0942a Mon Sep 17 00:00:00 2001
From: Joseph Weston <joseph.weston08@gmail.com>
Date: Thu, 11 May 2017 14:48:44 +0200
Subject: [PATCH] kpm: refactor the SpectralDensity API
It is now impossible to reduce the number of random vectors or moments.
It is also no longer possible to set the number of energy points used by
the method -- this can be supplied by providing a sequence of energies
when calling the KPM object. We replace the 'increase_energy_resolution'
and 'increase_accuracy' methods by 'add_moments' and 'add_vectors'.
The tests have also been updated to reflect these changes.
---
kwant/kpm.py | 222 +++++++++++++++++++---------------------
kwant/tests/test_kpm.py | 106 +++++++------------
2 files changed, 143 insertions(+), 185 deletions(-)
diff --git a/kwant/kpm.py b/kwant/kpm.py
index 2f4f8ac6..f55b532b 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -7,10 +7,7 @@
# the file AUTHORS.rst at the top-level directory of this distribution and at
# http://kwant-project.org/authors.
-__all__ = ['SpectralDensity']
-
-import warnings
-
+import math
import numpy as np
import scipy
import scipy.sparse.linalg as sla
@@ -20,6 +17,9 @@ from . import system
from ._common import ensure_rng
from .operator import _LocalOperator
+__all__ = ['SpectralDensity']
+
+
class SpectralDensity:
"""Calculate the spectral density of an operator.
@@ -50,13 +50,14 @@ class SpectralDensity:
``scipy.sparse.matrices``, the densities will be scalars.
If no operator is provided, the density of states is calculated
with a faster algorithm.
- num_vectors : integer, default: 10
+ num_vectors : positive int, default: 10
Number of random vectors for the KPM method.
- num_moments : integer, default: 100
- Number of moments, order of the KPM expansion.
- num_energies : integer, optional
- Number of points where the spectral density will be evaluated.
- If not provided, ``2*num_moments`` will be used.
+ num_moments : positive int, default: 100
+ Number of moments, order of the KPM expansion. Mutually exclusive
+ with 'energy_resolution'.
+ energy_resolution : positive float, optional
+ The resolution in energy of the KPM approximation to the spectral
+ density. Mutually exclusive with 'num_moments'.
vector_factory : function, optional
The user defined function ``f(n)`` generates random vectors of
length ``n`` that will be used in the algorithm.
@@ -66,7 +67,7 @@ class SpectralDensity:
bounds : pair of floats, optional
Lower and upper bounds for the eigenvalue spectrum of the system.
If not provided, they are computed.
- eps : float, default: 0.05
+ eps : positive float, default: 0.05
Parameter to ensure that the rescaled spectrum lies in the
interval ``(-1, 1)``; required for stability.
rng : seed, or random number generator, optional
@@ -126,15 +127,33 @@ class SpectralDensity:
Attributes
----------
energies : array of floats
- Array of sampling points with length ``num_energies`` in
+ Array of sampling points with length ``2 * num_moments`` in
the range of the spectrum.
densities : array of floats
Spectral density of the ``operator`` evaluated at the energies.
"""
def __init__(self, hamiltonian, params=None, operator=None,
- num_vectors=10, num_moments=100, num_energies=None,
+ num_vectors=10, num_moments=None, energy_resolution=None,
vector_factory=None, bounds=None, eps=0.05, rng=None):
+
+ if num_moments and energy_resolution:
+ raise TypeError("either 'num_moments' or 'energy_resolution' "
+ "must be provided.")
+
+ if energy_resolution:
+ num_moments = math.ceil((1.6 * self._a) / energy_resolution)
+ elif num_moments is None:
+ num_moments = 100
+
+ must_be_positive_int = ['num_vectors', 'num_moments']
+ for var in must_be_positive_int:
+ val = locals()[var]
+ if val <= 0 or val != int(val):
+ raise ValueError('{} must be a positive integer'.format(var))
+ if eps <= 0:
+ raise ValueError('eps must be positive')
+
rng = ensure_rng(rng)
# self.eps ensures that the rescaled Hamiltonian has a
# spectrum strictly in the interval (-1,1).
@@ -164,16 +183,6 @@ class SpectralDensity:
raise ValueError('Parameter `operator` has no `.dot` '
'attribute and is not callable.')
- self.num_moments = num_moments
- # Default number of sampling points
- if num_energies is None:
- self.num_energies = 2 * self.num_moments
- elif num_energies < self.num_moments:
- raise ValueError('The number of sampling points should be larger '
- 'than the number of moments.')
- else:
- self.num_energies = num_energies
-
self._vector_factory = vector_factory or \
(lambda n: np.exp(2j * np.pi * rng.random_sample(n)))
# store this vector for reproducibility
@@ -192,25 +201,15 @@ class SpectralDensity:
self._last_two_alphas = [0.] * num_vectors
self._moments_list = [0.] * num_vectors
- self.num_vectors = 0 # new random vectors will be used
+ self.num_moments = num_moments
+ self.num_vectors = 0 # new random vectors will be used
self._update_moments_list(self.num_moments, num_vectors)
- #update the number of random vectors
self.num_vectors = num_vectors
- # sum moments of all random vectors
- moments = np.sum(np.asarray(self._moments_list).real, axis=0)
- # divide by the number of random vectors
- moments = moments / self.num_vectors
- # divide by the length of the vectors to normalize
- moments = moments / self.hamiltonian.shape[0]
- # divide by scale factor to reflect the integral rescaling
- moments = moments / self._a
-
- # obtain energies, densities, and gammas of Eq. 83
+ moments = self._moments()
xk_rescaled, rho, self._gammas = _calc_fft_moments(
- moments, self.num_energies)
- # energies to normal scale
- self.energies = xk_rescaled*self._a + self._b
+ moments, 2 * self.num_moments)
+ self.energies = xk_rescaled * self._a + self._b
self.densities = rho
def __call__(self, energy=None):
@@ -239,15 +238,7 @@ class SpectralDensity:
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).real, axis=0)
- # divide by the number of random vectors
- moments = moments / self.num_vectors
- # divide by the length of the vectors to normalize
- moments = moments / self.hamiltonian.shape[0]
- # divide by scale factor to reflect the integral rescaling
- moments = moments / self._a
+ moments = self._moments()
m = np.arange(self.num_moments)
@@ -275,7 +266,7 @@ class SpectralDensity:
# 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_energies
+ factor = self._a / (2 * self.num_moments)
if distribution_function is None:
rho = self._gammas
else:
@@ -285,67 +276,79 @@ class SpectralDensity:
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``.
+ def add_moments(self, num_moments=None, *, energy_resolution=None):
+ """Increase the number of Chebyshev moments.
- The energy resolution is increased by increasing the number of
- sampling points. By default the number of moments is also
- increased.
+ Parameters
+ ----------
+ num_moments: positive int
+ The number of Chebyshev moments to add. Mutually
+ exclusive with 'energy_resolution'.
+ energy_resolution: positive float, optional
+ Features wider than this resolution are visible
+ in the spectral density. Mutually exclusive with
+ 'num_moments'.
"""
- resolution = 2 * self._a / (self.num_energies / 1.6)
- if tol > resolution:
- warnings.warn('Energy resolution is already smaller than tol.')
- else:
- num_energies = np.ceil((1.6 * 2*self._a) / tol)
- num_energies = num_energies.astype(int)
- if increase_num_moments:
- num_moments = num_energies // 2
- self.increase_accuracy(num_moments=num_moments,
- num_energies=num_energies)
- else:
- self.increase_accuracy(num_energies=num_energies)
+ if not ((num_moments is None) ^ (energy_resolution is None)):
+ raise TypeError("either 'num_moments' or 'energy_resolution' "
+ "must be provided.")
+
+ if energy_resolution:
+ if energy_resolution <= 0:
+ raise ValueError("'energy_resolution' must be positive"
+ .format(energy_resolution))
+ # factor of 1.6 comes from the fact that we use the
+ # Jackson kernel when calculating the FFT, which has
+ # maximal slope π/2. Rounding to 1.6 ensures that the
+ # energy resolution is sufficient.
+ present_resolution = self._a * 1.6 / self.num_moments
+ if present_resolution < energy_resolution:
+ raise ValueError('Energy resolution is already smaller '
+ 'than the requested resolution')
+ num_moments = math.ceil((1.6 * self._a) / energy_resolution)
+
+ if (num_moments is None or num_moments <= 0
+ or num_moments != int(num_moments)):
+ raise ValueError("'num_moments' must be a positive integer")
+
+ self._update_moments_list(self.num_moments + num_moments,
+ self.num_vectors)
+ self.num_moments += num_moments
+
+ # recalculate quantities derived from the moments
+ moments = self._moments()
+ xk_rescaled, rho, self._gammas = _calc_fft_moments(
+ moments, 2 * self.num_moments)
+ self.energies = xk_rescaled * self._a + self._b
+ self.densities = rho
- def increase_accuracy(self, num_moments=None, num_vectors=None,
- num_energies=None):
- """Increase the number of moments, random vectors, or sampling
- points.
+ def add_vectors(self, num_vectors):
+ """Increase the number of random vectors.
Parameters
----------
- num_moments, num_vectors, num_energies : integer
- Number of Chebyshev moments, random vectors used for
- sampling, and sampling points for the calculation of
- densities.
+ num_vectors: positive int
+ The number of random vectors to add.
"""
- # 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_vectors = num_vectors or self.num_vectors
- new_rand_vect = num_vectors - self.num_vectors
- for r in range(new_rand_vect):
+ if num_vectors <= 0 or num_vectors != int(num_vectors):
+ raise ValueError("'num_vectors' must be a positive integer")
+ for r in range(num_vectors):
self._rand_vect_list.append(
self._vector_factory(self.hamiltonian.shape[0]))
- self._moments_list.extend([0.] * new_rand_vect)
- self._last_two_alphas.extend([0.] * new_rand_vect)
- self._update_moments_list(self.num_moments, num_vectors)
- self.num_vectors = num_vectors
-
- num_moments = num_moments or self.num_moments
- self._update_moments_list(num_moments, self.num_vectors)
- self.num_moments = num_moments
-
- ## cut random vectors and moments when necessary
- self._moments_list[:] = [m[:num_moments]
- for m in self._moments_list[:num_vectors]]
-
- num_energies = num_energies or self.num_energies
- if num_energies < self.num_moments:
- raise ValueError(
- 'The number of sampling points should be larger than the '
- 'number of moments.')
- self.num_energies = num_energies
+ self._moments_list.extend([0.] * num_vectors)
+ self._last_two_alphas.extend([0.] * num_vectors)
+ self._update_moments_list(self.num_moments,
+ self.num_vectors + num_vectors)
+ self.num_vectors += num_vectors
+
+ # recalculate quantities derived from the moments
+ moments = self._moments()
+ xk_rescaled, rho, self._gammas = _calc_fft_moments(
+ moments, 2 * self.num_moments)
+ self.energies = xk_rescaled * self._a + self._b
+ self.densities = rho
+ def _moments(self):
# sum moments of all random vectors
moments = np.sum(np.asarray(self._moments_list).real, axis=0)
# divide by the number of random vectors
@@ -353,14 +356,7 @@ class SpectralDensity:
# divide by the length of the vectors to normalize
moments = moments / self.hamiltonian.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_energies)
-
- # Sampling points to normal scale
- self.energies = xk_rescaled * self._a + self._b
- self.densities = rho
+ return moments / self._a
def _update_moments_list(self, n_moments, n_rand):
"""Calculate the Chebyshev moments of an operator's spectral
@@ -385,8 +381,7 @@ class SpectralDensity:
r_start = self.num_vectors
new_rand_vect = n_rand - self.num_vectors
else:
- warnings.warn('Decreasing number of random vectors')
- return
+ raise ValueError('Cannot decrease number of random vectors')
if n_moments == self.num_moments:
m_start = 2
@@ -398,12 +393,11 @@ class SpectralDensity:
new_moments = n_moments - self.num_moments
m_start = self.num_moments
if new_moments < 0:
- warnings.warn('Decreasing the number of moments.')
- return
+ raise ValueError('Cannot decrease number of moments')
- assert new_rand_vect == 0,\
- ("Only 'num_moments' *or* 'num_vectors' may be updated "
- "at a time.")
+ if new_rand_vect != 0:
+ raise ValueError("Only 'num_moments' *or* 'num_vectors' "
+ "may be updated at a time.")
for r in range(r_start, n_rand):
alpha_zero = self._rand_vect_list[r]
diff --git a/kwant/tests/test_kpm.py b/kwant/tests/test_kpm.py
index 94b411d5..5ea3dbf7 100644
--- a/kwant/tests/test_kpm.py
+++ b/kwant/tests/test_kpm.py
@@ -25,7 +25,6 @@ dim = 20
p = SimpleNamespace(
num_moments=200,
num_vectors=5,
- num_energies=400
)
ham = kwant.rmt.gaussian(dim)
@@ -51,7 +50,6 @@ def make_spectrum(ham, p, operator=None, vector_factory=None, rng=None, params=N
vector_factory=vector_factory,
num_moments=p.num_moments,
num_vectors=p.num_vectors,
- num_energies=p.num_energies,
rng=rng,
params=params
)
@@ -172,42 +170,6 @@ def test_api_operator():
operator=not_an_operator)
-def test_api_lower_num_moments():
- spectrum = SpectralDensity(kwant.rmt.gaussian(dim))
- num_moments = spectrum.num_moments - 1
- with pytest.warns(UserWarning):
- 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_energies=num_samplint_points)
-
-
-def test_api_warning_less_sampling_points():
- precise = copy(p)
- precise.num_energies = 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_vectors():
- spectrum = SpectralDensity(kwant.rmt.gaussian(dim))
- num_vectors = spectrum.num_vectors - 1
- with pytest.warns(UserWarning):
- spectrum.increase_accuracy(num_vectors=num_vectors)
-
-
def test_api_single_eigenvalue_error():
with pytest.raises(ValueError):
SpectralDensity(np.identity(dim, dtype=complex))
@@ -248,6 +210,7 @@ def test_bounds():
# 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.
@@ -357,7 +320,7 @@ def test_increase_num_moments():
precise.num_moments = 2 * p.num_moments
spectrum_raise = make_spectrum(ham, precise, rng=1)
- spectrum.increase_accuracy(num_moments=precise.num_moments)
+ spectrum.add_moments(p.num_moments)
# Check bit for bit equality
assert np.all(np.array(spectrum_raise._moments_list) ==
@@ -373,13 +336,11 @@ def test_increase_num_moments():
# test when increasing num_moments by 1 from an even to an odd number
spectrum_even = make_spectrum(ham, p, rng=1)
- spectrum_even.increase_accuracy(num_moments=p.num_moments + 1)
+ spectrum_even.add_moments(num_moments=1)
assert np.all(np.array(spectrum_even._moments_list) ==
np.array(spectrum_odd._moments_list))
-
-
# ### increase num_moments with an operator (different algorithm)
def test_increase_num_moments_op():
@@ -389,8 +350,8 @@ def test_increase_num_moments_op():
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)
+ sp1.add_moments(num_moments=sp1.num_moments)
+ sp2.add_moments(num_moments=sp2.num_moments)
# different algorithms are used so these arrays are equal up to TOL_SP
assert_allclose_sp(sp1.densities, sp2.densities)
@@ -404,24 +365,41 @@ def test_increase_num_vectors():
spectrum_raise = make_spectrum(ham, precise, rng=1)
spectrum = make_spectrum(ham, p, rng=1)
- spectrum.increase_accuracy(num_vectors=precise.num_vectors)
+ spectrum.add_vectors(num_vectors=p.num_vectors)
# Check bit for bit equality
assert np.all(spectrum_raise.densities == spectrum.densities)
-# ### increase num_energies
+def test_invalid_input():
-def test_increase_num_energies():
- precise = copy(p)
- precise.sampling_points = 2 * p.num_energies
+ with pytest.raises(TypeError):
+ SpectralDensity(ham, num_moments=10, energy_resolution=0.1)
- spectrum_raise = make_spectrum(ham, precise, rng=1)
- spectrum = make_spectrum(ham, p, rng=1)
- spectrum.increase_accuracy(num_moments=precise.num_moments)
+ s = SpectralDensity(ham)
+ with pytest.raises(TypeError):
+ s.add_moments(num_moments=20, energy_resolution=0.1)
+
+ for v in (-20, 0, 10.5):
+ for p in ('num_vectors', 'num_moments'):
+ with pytest.raises(ValueError):
+ SpectralDensity(ham, **{p: v})
+ s = SpectralDensity(ham)
+ with pytest.raises(ValueError):
+ s.add_moments(num_moments=v)
+ with pytest.raises(ValueError):
+ s.add_vectors(num_vectors=v)
+
+ for v in (-1, -0.5, 0):
+ with pytest.raises(ValueError):
+ SpectralDensity(ham, eps=v)
+ s = SpectralDensity(ham)
+ with pytest.raises(ValueError):
+ s.add_moments(energy_resolution=v)
+
+ with pytest.raises(ValueError):
+ s.add_moments(energy_resolution=10)
- # Check bit for bit equality
- assert np.all(spectrum_raise.densities == spectrum.densities)
# ### Convergence for higher number of moments ###
@@ -434,7 +412,6 @@ def test_check_convergence_decreasing_values():
precise = SimpleNamespace(
num_moments=10*dim + 100*i*dim,
num_vectors=dim//2,
- num_energies=None
)
results = []
np.random.seed(1)
@@ -480,7 +457,6 @@ def test_convergence_custom_vector_factory():
precise = SimpleNamespace(
num_moments=10*dim + 100*i*dim,
num_vectors=dim//2,
- num_energies=None
)
results = []
iterations = 3
@@ -537,28 +513,16 @@ def test_average():
def test_increase_energy_resolution():
spectrum, filter_index = make_spectrum_and_peaks(ham, p)
- old_sampling_points = spectrum.num_energies
+ old_sampling_points = 2 * spectrum.num_moments
tol = np.max(np.abs(np.diff(spectrum.energies))) / 2
- spectrum.increase_energy_resolution(tol=tol)
- new_sampling_points = spectrum.num_energies
+ spectrum.add_moments(energy_resolution=tol)
+ new_sampling_points = 2 * spectrum.num_moments
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_energies
- 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_energies
-
- assert old_sampling_points < new_sampling_points
- assert np.max(np.abs(np.diff(spectrum.energies))) < tol
-
# ### check _rescale
--
GitLab