From bb8edcd5a220b74a82645957f2e774692e2e5caf Mon Sep 17 00:00:00 2001
From: Joseph Weston <joseph@weston.cloud>
Date: Mon, 9 Sep 2019 13:41:19 +0200
Subject: [PATCH] update tests to use new finalized Builder interface

---
 kwant/tests/test_builder.py | 194 +++++++++++++++++++++++++-----------
 kwant/tests/test_system.py  |  52 ++++++----
 2 files changed, 167 insertions(+), 79 deletions(-)

diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index f41652a8..f9832e38 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -1,4 +1,4 @@
-# Copyright 2011-2018 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
@@ -290,12 +290,28 @@ def random_hopping_integral(rng):
 
 
 def check_onsite(fsyst, sites, subset=False, check_values=True):
+    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
+
+    if vectorized:
+        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
+
     freq = {}
     for node in range(fsyst.graph.num_nodes):
         site = fsyst.sites[node].tag
         freq[site] = freq.get(site, 0) + 1
         if check_values and site in sites:
-            assert fsyst.onsites[node][0] is sites[site]
+            if vectorized:
+                term = fsyst._onsite_term_by_site_id[node]
+                value = fsyst._term_values[term]
+                if callable(value):
+                    assert value is sites[site]
+                else:
+                    (w, _), (off, _) = fsyst.subgraphs[fsyst.terms[term].subgraph]
+                    node_off = node - site_offsets[w]
+                    selector = np.searchsorted(off, node_off)
+                    assert np.allclose(value[selector], sites[site])
+            else:
+                assert fsyst.onsites[node][0] is sites[site]
     if not subset:
         # Check that all sites of `fsyst` are in `sites`.
         for site in freq.keys():
@@ -306,24 +322,50 @@ def check_onsite(fsyst, sites, subset=False, check_values=True):
 
 
 def check_hoppings(fsyst, hops):
+    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
+
+    if vectorized:
+        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
+
     assert fsyst.graph.num_edges == 2 * len(hops)
     for edge_id, edge in enumerate(fsyst.graph):
-        tail, head = edge
-        tail = fsyst.sites[tail].tag
-        head = fsyst.sites[head].tag
-        value = fsyst.hoppings[edge_id][0]
-        if value is builder.Other:
-            assert (head, tail) in hops
+        i, j = edge
+        tail = fsyst.sites[i].tag
+        head = fsyst.sites[j].tag
+
+        if vectorized:
+            term = fsyst._hopping_term_by_edge_id[edge_id]
+            if term < 0:  # Hermitian conjugate
+                assert (head, tail) in hops
+            else:
+                value = fsyst._term_values[term]
+                assert (tail, head) in hops
+                if callable(value):
+                    assert value is hops[tail, head]
+                else:
+                    dtype = np.dtype([('f0', int), ('f1', int)])
+                    subgraph = fsyst.terms[term].subgraph
+                    (to_w, from_w), hoppings = fsyst.subgraphs[subgraph]
+                    hop = (i - site_offsets[to_w], j - site_offsets[from_w])
+                    hop = np.array(hop, dtype=dtype)
+                    hoppings = hoppings.transpose().view(dtype).reshape(-1)
+                    selector = np.recarray.searchsorted(hoppings, hop)
+                    assert np.allclose(value[selector], hops[tail, head])
         else:
-            assert (tail, head) in hops
-            assert value is hops[tail, head]
+            value = fsyst.hoppings[edge_id][0]
+            if value is builder.Other:
+                assert (head, tail) in hops
+            else:
+                assert (tail, head) in hops
+                assert value is hops[tail, head]
 
 def check_id_by_site(fsyst):
     for i, site in enumerate(fsyst.sites):
         assert fsyst.id_by_site[site] == i
 
 
-def test_finalization():
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_finalization(vectorize):
     """Test the finalization of finite and infinite systems.
 
     In order to exactly verify the finalization, low-level features of the
@@ -377,7 +419,7 @@ def test_finalization():
     neighbors = sorted(neighbors)
 
     # Build scattering region from blueprint and test it.
-    syst = builder.Builder()
+    syst = builder.Builder(vectorize=vectorize)
     fam = kwant.lattice.general(ta.identity(2), norbs=1)
     for site, value in sr_sites.items():
         syst[fam(*site)] = value
@@ -388,7 +430,7 @@ def test_finalization():
     check_onsite(fsyst, sr_sites)
     check_hoppings(fsyst, sr_hops)
     # check that sites are sorted
-    assert fsyst.sites == tuple(sorted(fam(*site) for site in sr_sites))
+    assert tuple(fsyst.sites) == tuple(sorted(fam(*site) for site in sr_sites))
 
     # Build lead from blueprint and test it.
     lead = builder.Builder(kwant.TranslationalSymmetry((size, 0)))
@@ -438,7 +480,7 @@ def test_finalization():
     # site ordering independent of whether interface was specified
     flead_order = builder.InfiniteSystem(lead, [system.Site(fam, n)
                                                 for n in neighbors])
-    assert flead.sites == flead_order.sites
+    assert tuple(flead.sites) == tuple(flead_order.sites)
 
 
     syst.leads[-1] = builder.BuilderLead(
@@ -457,47 +499,62 @@ def test_finalization():
     raises(ValueError, lead.finalized)
 
 
-def test_site_ranges():
+def _make_system_from_sites(sites, vectorize):
+    syst = builder.Builder(vectorize=vectorize)
+    for s in sites:
+        norbs = s.family.norbs
+        if norbs:
+            syst[s] = np.eye(s.family.norbs)
+        else:
+            syst[s] = None
+    return syst.finalized()
+
+
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_site_ranges(vectorize):
     lat1a = kwant.lattice.chain(norbs=1, name='a')
     lat1b = kwant.lattice.chain(norbs=1, name='b')
     lat2 = kwant.lattice.chain(norbs=2)
-    site_ranges = builder._site_ranges
 
     # simple case -- single site family
     for lat in (lat1a, lat2):
         sites = list(map(lat, range(10)))
-        ranges = site_ranges(sites)
+        syst = _make_system_from_sites(sites, vectorize)
+        ranges = syst.site_ranges
         expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
-        assert ranges == expected
+        assert np.array_equal(ranges, expected)
 
     # pair of site families
-    sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)),
-                     map(lat1a, range(4)))
-    expected = [(0, 1, 0), (4, 1, 4), (10, 1, 10), (14, 0, 14)]
-    assert expected == site_ranges(tuple(sites))
+    sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)))
+    syst = _make_system_from_sites(sites, vectorize)
+    expected = [(0, 1, 0), (4, 1, 4), (10, 0, 10)]
+    assert np.array_equal(expected, syst.site_ranges)
     sites = it.chain(map(lat2, range(4)), map(lat1a, range(6)),
                      map(lat1b, range(4)))
+    syst = _make_system_from_sites(sites, vectorize)
     expected = [(0, 2, 0), (4, 1, 4*2), (10, 1, 4*2+6), (14, 0, 4*2+10)]
-    assert expected == site_ranges(tuple(sites))
+    assert np.array_equal(expected, syst.site_ranges)
 
     # test with an actual builder
     for lat in (lat1a, lat2):
         sites = list(map(lat, range(10)))
-        syst = kwant.Builder()
+        syst = kwant.Builder(vectorize=vectorize)
         syst[sites] = np.eye(lat.norbs)
         ranges = syst.finalized().site_ranges
         expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
-        assert ranges == expected
-        # poison system with a single site with no norbs defined.
-        # Also catch the deprecation warning.
-        with warnings.catch_warnings():
-            warnings.simplefilter("ignore")
-            syst[kwant.lattice.chain(norbs=None)(0)] = 1
-        ranges = syst.finalized().site_ranges
-        assert ranges == None
-
-
-def test_hamiltonian_evaluation():
+        assert np.array_equal(ranges, expected)
+        if not vectorize:  # vectorized systems *must* have all norbs
+            # poison system with a single site with no norbs defined.
+            # Also catch the deprecation warning.
+            with warnings.catch_warnings():
+                warnings.simplefilter("ignore")
+                syst[kwant.lattice.chain(norbs=None)(0)] = 1
+            ranges = syst.finalized().site_ranges
+            assert ranges is None
+
+
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_hamiltonian_evaluation(vectorize):
     def f_onsite(site):
         return site.tag[0]
 
@@ -505,14 +562,26 @@ def test_hamiltonian_evaluation():
         a, b = a.tag, b.tag
         return complex(a[0] + b[0], a[1] - b[1])
 
+    def f_onsite_vectorized(sites):
+        return sites.tags[:, 0]
+
+    def f_hopping_vectorized(a, b):
+        a, b = a.tags, b.tags
+        return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])
+
     tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
     edges = [(0, 1), (0, 2), (0, 3), (1, 2)]
 
-    syst = builder.Builder()
+    syst = builder.Builder(vectorize=vectorize)
     fam = builder.SimpleSiteFamily(norbs=1)
     sites = [fam(*tag) for tag in tags]
-    syst[(fam(*tag) for tag in tags)] = f_onsite
-    syst[((fam(*tags[i]), fam(*tags[j])) for (i, j) in edges)] = f_hopping
+    hoppings = [(sites[i], sites[j]) for i, j in edges]
+    if vectorize:
+        syst[sites] = f_onsite_vectorized
+        syst[hoppings] = f_hopping_vectorized
+    else:
+        syst[sites] = f_onsite
+        syst[hoppings] = f_hopping
     fsyst = syst.finalized()
 
     assert fsyst.graph.num_nodes == len(tags)
@@ -521,12 +590,12 @@ def test_hamiltonian_evaluation():
     for i in range(len(tags)):
         site = fsyst.sites[i]
         assert site in sites
-        assert fsyst.hamiltonian(i, i) == syst[site](site)
+        assert fsyst.hamiltonian(i, i) == f_onsite(site)
 
     for t, h in fsyst.graph:
         tsite = fsyst.sites[t]
         hsite = fsyst.sites[h]
-        assert fsyst.hamiltonian(t, h) == syst[tsite, hsite](tsite, hsite)
+        assert fsyst.hamiltonian(t, h) == f_hopping(tsite, hsite)
 
     # test when user-function raises errors
     def onsite_raises(site):
@@ -555,7 +624,7 @@ def test_hamiltonian_evaluation():
     syst[new_hop[0]] = onsite_raises
     syst[new_hop] = hopping_raises
     fsyst = syst.finalized()
-    hop = tuple(map(fsyst.sites.index, new_hop))
+    hop = tuple(map(fsyst.id_by_site.__getitem__, new_hop))
     test_raising(fsyst, hop)
 
     # test with infinite system
@@ -563,7 +632,7 @@ def test_hamiltonian_evaluation():
     for k, v in it.chain(syst.site_value_pairs(), syst.hopping_value_pairs()):
         inf_syst[k] = v
     inf_fsyst = inf_syst.finalized()
-    hop = tuple(map(inf_fsyst.sites.index, new_hop))
+    hop = tuple(map(inf_fsyst.id_by_site.__getitem__, new_hop))
     test_raising(inf_fsyst, hop)
 
 
@@ -1204,15 +1273,22 @@ def test_discrete_symmetries():
 # We need to keep testing 'args', but we don't want to see
 # all the deprecation warnings in the test logs
 @pytest.mark.filterwarnings("ignore:.*'args' parameter")
-def test_argument_passing():
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_argument_passing(vectorize):
     chain = kwant.lattice.chain(norbs=1)
 
     # Test for passing parameters to hamiltonian matrix elements
     def onsite(site, p1, p2):
-        return p1 + p2
+        if vectorize:
+            return (p1 + p2) * np.ones(len(site))
+        else:
+            return p1 + p2
 
     def hopping(site1, site2, p1, p2):
-        return p1 - p2
+        if vectorize:
+            return (p1 - p2) * np.ones(len(site1))
+        else:
+            return p1 - p2
 
     def gen_fill_syst(onsite, hopping, syst):
         syst[(chain(i) for i in range(3))] = onsite
@@ -1221,8 +1297,9 @@ def test_argument_passing():
 
     fill_syst = ft.partial(gen_fill_syst, onsite, hopping)
 
-    syst = fill_syst(kwant.Builder())
-    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
+    syst = fill_syst(kwant.Builder(vectorize=vectorize))
+    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)),
+                                       vectorize=vectorize))
 
     tests = (
         syst.hamiltonian_submatrix,
@@ -1250,16 +1327,12 @@ def test_argument_passing():
 
     # test that passing parameters without default values works, and that
     # passing parameters with default values fails
-    def onsite(site, p1, p2):
-        return p1 + p2
-
-    def hopping(site, site2, p1, p2):
-        return p1 - p2
 
     fill_syst = ft.partial(gen_fill_syst, onsite, hopping)
 
-    syst = fill_syst(kwant.Builder())
-    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
+    syst = fill_syst(kwant.Builder(vectorize=vectorize))
+    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)),
+                                       vectorize=vectorize))
 
     tests = (
         syst.hamiltonian_submatrix,
@@ -1275,12 +1348,18 @@ def test_argument_passing():
 
     # Some common, some different args for value functions
     def onsite2(site, a, b):
-        return site.pos[0] + a + b
+        if vectorize:
+            return site.positions()[:, 0] + a + b
+        else:
+            return site.pos[0] + a + b
 
     def hopping2(site1, site2, a, c, b):
-        return a + b + c
+        if vectorize:
+            return (a + b + c) * np.ones(len(site1))
+        else:
+            return a + b + c
 
-    syst = kwant.Builder()
+    syst = kwant.Builder(vectorize=vectorize)
     syst[(chain(i) for i in range(3))] = onsite2
     syst[((chain(i), chain(i + 1)) for i in range(2))] = hopping2
     fsyst = syst.finalized()
@@ -1391,6 +1470,7 @@ def test_subs():
     lead = lead.finalized()
     assert lead.parameters == {'lead_a', 'lead_b', 'lead_c'}
 
+
 def test_attach_stores_padding():
     lat = kwant.lattice.chain(norbs=1)
     syst = kwant.Builder()
diff --git a/kwant/tests/test_system.py b/kwant/tests/test_system.py
index 9bba152b..3a403477 100644
--- a/kwant/tests/test_system.py
+++ b/kwant/tests/test_system.py
@@ -1,4 +1,4 @@
-# Copyright 2011-2016 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
@@ -8,6 +8,7 @@
 
 import pickle
 import copy
+import pytest
 from pytest import raises
 import numpy as np
 from scipy import sparse
@@ -15,9 +16,11 @@ import kwant
 from kwant._common import ensure_rng
 
 
-def test_hamiltonian_submatrix():
-    syst = kwant.Builder()
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_hamiltonian_submatrix(vectorize):
+    syst = kwant.Builder(vectorize=vectorize)
     chain = kwant.lattice.chain(norbs=1)
+    chain2 = kwant.lattice.chain(norbs=2)
     for i in range(3):
         syst[chain(i)] = 0.5 * i
     for i in range(2):
@@ -27,7 +30,11 @@ def test_hamiltonian_submatrix():
     mat = syst2.hamiltonian_submatrix()
     assert mat.shape == (3, 3)
     # Sorting is required due to unknown compression order of builder.
-    perm = np.argsort([os[0] for os in syst2.onsites])
+    if vectorize:
+        _, (site_offsets, _) = syst2.subgraphs[0]
+    else:
+        site_offsets = [os[0] for os in syst2.onsites]
+    perm = np.argsort(site_offsets)
     mat_should_be = np.array([[0, 1j, 0], [-1j, 0.5, 2j], [0, -2j, 1]])
 
     mat = mat[perm, :]
@@ -41,19 +48,12 @@ def test_hamiltonian_submatrix():
     mat = mat[:, perm]
     np.testing.assert_array_equal(mat, mat_should_be)
 
-    mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]])
-    np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
-
-    mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]], sparse=True)
-    mat = mat.toarray()
-    np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
-
     # Test for correct treatment of matrix input.
-    syst = kwant.Builder()
-    syst[chain(0)] = np.array([[0, 1j], [-1j, 0]])
+    syst = kwant.Builder(vectorize=vectorize)
+    syst[chain2(0)] = np.array([[0, 1j], [-1j, 0]])
     syst[chain(1)] = np.array([[1]])
     syst[chain(2)] = np.array([[2]])
-    syst[chain(1), chain(0)] = np.array([[1, 2j]])
+    syst[chain(1), chain2(0)] = np.array([[1, 2j]])
     syst[chain(2), chain(1)] = np.array([[3j]])
     syst2 = syst.finalized()
     mat_dense = syst2.hamiltonian_submatrix()
@@ -62,9 +62,10 @@ def test_hamiltonian_submatrix():
 
     # Test precalculation of modes.
     rng = ensure_rng(5)
-    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
-    lead[chain(0)] = np.zeros((2, 2))
-    lead[chain(0), chain(1)] = rng.randn(2, 2)
+    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)),
+                         vectorize=vectorize)
+    lead[chain2(0)] = np.zeros((2, 2))
+    lead[chain2(0), chain2(1)] = rng.randn(2, 2)
     syst.attach_lead(lead)
     syst2 = syst.finalized()
     smatrix = kwant.smatrix(syst2, .1).data
@@ -74,19 +75,26 @@ def test_hamiltonian_submatrix():
     raises(ValueError, kwant.solvers.default.greens_function, syst3, 0.2)
 
     # Test for shape errors.
-    syst[chain(0), chain(2)] = np.array([[1, 2]])
+    syst[chain2(0), chain(2)] = np.array([[1, 2]])
     syst2 = syst.finalized()
     raises(ValueError, syst2.hamiltonian_submatrix)
     raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
-    syst[chain(0), chain(2)] = 1
+    syst[chain2(0), chain(2)] = 1
     syst2 = syst.finalized()
     raises(ValueError, syst2.hamiltonian_submatrix)
     raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
+    if vectorize:  # non-vectorized systems don't check this at finalization
+        # Add another hopping of the same type but with a different
+        # (and still incompatible) shape.
+        syst[chain2(0), chain(1)] = np.array([[1, 2]])
+        raises(ValueError, syst.finalized)
 
 
-def test_pickling():
-    syst = kwant.Builder()
-    lead = kwant.Builder(symmetry=kwant.TranslationalSymmetry([1.]))
+@pytest.mark.parametrize("vectorize", [False, True])
+def test_pickling(vectorize):
+    syst = kwant.Builder(vectorize=vectorize)
+    lead = kwant.Builder(symmetry=kwant.TranslationalSymmetry([1.]),
+                         vectorize=vectorize)
     lat = kwant.lattice.chain(norbs=1)
     syst[lat(0)] = syst[lat(1)] = 0
     syst[lat(0), lat(1)] = 1
-- 
GitLab