From d5eb75ca2d486e6dc779f9edce6050f899eb16e1 Mon Sep 17 00:00:00 2001
From: Pablo Piskunow <pablo.perez.piskunow@gmail.com>
Date: Thu, 15 Jun 2017 09:54:19 +0200
Subject: [PATCH] fix bug in passing `energy_resolution` instead of
`num_moments`
When passing `energy_resolution`, the number of moments is set using
the yet undefined parameter `self._a`. Now the number of moments is
set after defining `self._a`.
---
kwant/kpm.py | 26 +++++++++++++-------------
kwant/tests/test_kpm.py | 18 ++++++++++++++++++
2 files changed, 31 insertions(+), 13 deletions(-)
diff --git a/kwant/kpm.py b/kwant/kpm.py
index 595cba04..2c4ee432 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -141,19 +141,6 @@ class SpectralDensity:
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).
@@ -195,6 +182,19 @@ class SpectralDensity:
bounds=bounds)
self.bounds = (self._b - self._a, self._b + self._a)
+ 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')
+
for r in range(num_vectors):
self._rand_vect_list.append(
self._vector_factory(self.hamiltonian.shape[0]))
diff --git a/kwant/tests/test_kpm.py b/kwant/tests/test_kpm.py
index 23529efc..c5b8806b 100644
--- a/kwant/tests/test_kpm.py
+++ b/kwant/tests/test_kpm.py
@@ -175,6 +175,24 @@ def test_api_single_eigenvalue_error():
SpectralDensity(np.identity(dim, dtype=complex))
+def test_energy_resolution():
+ """Check that energy resolution works and gives the same output as
+ setting the equivalent number of moments.
+ """
+ ham = kwant.rmt.gaussian(dim)
+ rng = 1
+ energy_resolution = 0.05
+
+ sp1 = SpectralDensity(ham, rng=rng, energy_resolution=energy_resolution)
+ # number of moments for this energy resolution
+ num_moments = sp1.num_moments
+ # use the same number of moments but not passing energy resolution
+ sp2 = SpectralDensity(ham, rng=rng, num_moments=num_moments)
+
+ # Check bit for bit equality
+ assert np.all(sp1.densities == sp2.densities)
+
+
def test_operator_none():
"""Check operator=None gives the same results as operator=np.identity(),
with the same random vectors.
--
GitLab