Skip to content
Snippets Groups Projects
Commit 3f675e33 authored by Joseph Weston's avatar Joseph Weston
Browse files

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.
parent 8e661bc9
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
kwant/kpm.py
+ 108
− 114
View file @ 3f675e33
......@@ -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]
......
kwant/tests/test_kpm.py
+ 35
− 71
View file @ 3f675e33
......@@ -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
......
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