Browse code

Merge branch 'feature/vectorization' into 'master'

vectorize systems so that value functions may be called once
on a whole vector of sites (for onsite value functions) or a pair of
vectors of sites (for hopping value functions). In this way we avoid
the Python overhead of calling the value functions many times. For systems
with vectorization enabled this can result in a ~100x speedup when evaluating
the Hamiltonian.

See merge request kwant/kwant!323
Closes kwant/kwant#302

Joseph Weston authored on 13/11/2019 16:33:41
Showing 23 changed files
... ...
@@ -42,7 +42,7 @@ description check `~kwant.plotter.plot` documentation.
42 42
 
43 43
 The behavior of `~kwant.plotter.plot` with low level systems has changed.
44 44
 Arguments of plot which are functions are given site numbers in place of
45
-`~kwant.builder.Site` objects when plotting a low level system.  This
45
+`~kwant.system.Site` objects when plotting a low level system.  This
46 46
 provides an easy way to make the appearance of lines and symbols depend on
47 47
 computation results.
48 48
 
... ...
@@ -39,7 +39,7 @@ now it suffices to write ::
39 39
     syst[lat.neighbors()] = t
40 40
 
41 41
 This is possible because `~kwant.builder.Builder` now accepts *functions* as
42
-keys in addition to `~kwant.builder.Site` objects and tuples of them
42
+keys in addition to `~kwant.system.Site` objects and tuples of them
43 43
 (hoppings).  These functions are expected to yield either sites or hoppings,
44 44
 when given a builder instance as the sole argument. The use of such keys is to
45 45
 implement sets of sites or hoppings that depend on what is already present in
... ...
@@ -3,6 +3,40 @@ What's new in Kwant 1.5
3 3
 
4 4
 This article explains the user-visible changes in Kwant 1.5.0.
5 5
 
6
+
7
+Value functions may now be vectorized
8
+-------------------------------------
9
+It is now possible to define value functions that act on whole
10
+arrays of sites at a time.
11
+::
12
+
13
+    def onsite(sites):
14
+        r = sites.positions()
15
+        return np.exp(-np.linalg.norm(r, axis=1)**2)
16
+
17
+    def hopping(to_sites, from_sites):
18
+        dr = to_sites.positions() - from_sites.positions()
19
+        return 1 / np.linalg.norm(dr, axis=1)**2
20
+
21
+
22
+    lat = kwant.lattice.square(norbs=1)
23
+    syst = kwant.Builder(vectorize=True)
24
+    syst[(lat(i, j) for i in range(4) for j in range(10)] = onsite
25
+    syst[lat.neighbors()] = hopping
26
+    syst = syst.finalized()
27
+
28
+    syst.hamiltonian_submatrix()
29
+
30
+Notice that when creating the Builder we provide the ``vectorize=True`` flag,
31
+and that the value functions now take `~kwant.system.SiteArray` objects
32
+(which have a ``positions`` method, rather than a ``pos`` property as for
33
+`~kwant.system.Site`). We can now operate on the array of site positions
34
+using ``numpy``. Using vectorized value functions in this way can produce an
35
+order of magnitude speedup when evaluating system Hamiltonians (though this
36
+speedup may be masked by other parts of your computation e.g. calculating
37
+the scattering matrix).
38
+
39
+
6 40
 Deprecation of leaving 'norbs' unset for site families
7 41
 ------------------------------------------------------
8 42
 When constructing site families (e.g. lattices) it is now deprecated to
... ...
@@ -9,7 +9,6 @@ Types
9 9
    :toctree: generated/
10 10
 
11 11
    Builder
12
-   Site
13 12
    HoppingKind
14 13
    SimpleSiteFamily
15 14
    BuilderLead
... ...
@@ -17,13 +16,14 @@ Types
17 16
    ModesLead
18 17
    FiniteSystem
19 18
    InfiniteSystem
19
+   FiniteVectorizedSystem
20
+   InfiniteVectorizedSystem
20 21
 
21 22
 Abstract base classes
22 23
 ---------------------
23 24
 .. autosummary::
24 25
    :toctree: generated/
25 26
 
26
-   SiteFamily
27 27
    Symmetry
28 28
    Lead
29 29
 
... ...
@@ -33,3 +33,11 @@ Functions
33 33
    :toctree: generated/
34 34
 
35 35
    add_peierls_phase
36
+
37
+Mixin Classes
38
+-------------
39
+.. autosummary::
40
+   :toctree: generated/
41
+
42
+   _FinalizedBuilderMixin
43
+   _VectorizedFinalizedBuilderMixin
... ...
@@ -19,10 +19,32 @@ In practice, very often Kwant systems are finalized
19 19
 `kwant.builder.FiniteSystem` or `kwant.builder.InfiniteSystem`) and offer some
20 20
 additional attributes.
21 21
 
22
+Sites
23
+-----
24
+.. autosummary::
25
+   :toctree: generated/
26
+
27
+   Site
28
+   SiteArray
29
+   SiteFamily
30
+
31
+Systems
32
+-------
22 33
 .. autosummary::
23 34
    :toctree: generated/
24 35
 
25 36
    System
37
+   VectorizedSystem
26 38
    InfiniteSystem
39
+   InfiniteVectorizedSystem
27 40
    FiniteSystem
41
+   FiniteVectorizedSystem
28 42
    PrecalculatedLead
43
+
44
+Mixin Classes
45
+-------------
46
+.. autosummary::
47
+   :toctree: generated/
48
+
49
+   FiniteSystemMixin
50
+   InfiniteSystemMixin
... ...
@@ -358,13 +358,13 @@ subbands that increases with energy.
358 358
      **sites** as indices. Sites themselves have a certain type, and
359 359
      belong to a **site family**. A site family is also used to convert
360 360
      something that represents a site (like a tuple) into a
361
-     proper `~kwant.builder.Site` object that can be used with
361
+     proper `~kwant.system.Site` object that can be used with
362 362
      `~kwant.builder.Builder`.
363 363
 
364 364
      In the above example, `lat` is the site family. ``lat(i, j)``
365 365
      then translates the description of a lattice site in terms of two
366 366
      integer indices (which is the natural way to do here) into
367
-     a proper `~kwant.builder.Site` object.
367
+     a proper `~kwant.system.Site` object.
368 368
 
369 369
      The concept of site families and sites allows `~kwant.builder.Builder`
370 370
      to mix arbitrary lattices and site families
... ...
@@ -438,7 +438,7 @@ the `~kwant.operator.Current`, instead of a constant matrix:
438 438
     # evaluate the operator
439 439
     m_current = J_m(psi, params=dict(r0=25, delta=10))
440 440
 
441
-The function must take a `~kwant.builder.Site` as its first parameter,
441
+The function must take a `~kwant.system.Site` as its first parameter,
442 442
 and may optionally take other parameters (i.e. it must have the same
443 443
 signature as a Hamiltonian onsite function), and must return the square
444 444
 matrix that defines the operator we wish to calculate.
... ...
@@ -529,8 +529,8 @@ We see that the probability current is conserved across the scattering region,
529 529
 but the z-projected spin current is not due to the fact that the Hamiltonian
530 530
 does not commute with :math:`σ_z` everywhere in the scattering region.
531 531
 
532
-.. note:: ``where`` can also be provided as a sequence of `~kwant.builder.Site`
533
-          or a sequence of hoppings (i.e. pairs of `~kwant.builder.Site`),
532
+.. note:: ``where`` can also be provided as a sequence of `~kwant.system.Site`
533
+          or a sequence of hoppings (i.e. pairs of `~kwant.system.Site`),
534 534
           rather than a function.
535 535
 
536 536
 
... ...
@@ -105,7 +105,7 @@ the start end end site of hopping as arguments:
105 105
     kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw);
106 106
 
107 107
 Note that since we are using an unfinalized Builder, a ``site`` is really an
108
-instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
108
+instance of `~kwant.system.Site`. With these adjustments we arrive at a plot
109 109
 that carries the same information, but is much easier to interpret.
110 110
 
111 111
 Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
... ...
@@ -175,7 +175,7 @@ this is interpreted as the radius of the inner circle).
175 175
 
176 176
 Finally, note that since we are dealing with a finalized system now, a site `i`
177 177
 is represented by an integer. In order to obtain the original
178
-`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
178
+`~kwant.system.Site`, ``syst.sites[i]`` can be used.
179 179
 
180 180
 The way how data is presented of course influences what features of the data
181 181
 are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
... ...
@@ -285,7 +285,7 @@ define the potential profile of a quantum well as:
285 285
         else:
286 286
             return 0
287 287
 
288
-This function takes two arguments: the first of type `~kwant.builder.Site`,
288
+This function takes two arguments: the first of type `~kwant.system.Site`,
289 289
 from which you can get the real-space coordinates using ``site.pos``, and the
290 290
 value of the potential as the second.  Note that in `potential` we can access
291 291
 variables `L` and `L_well` that are defined globally.
... ...
@@ -375,10 +375,10 @@ oscillatory transmission behavior through resonances in the quantum well.
375 375
 .. specialnote:: Technical details
376 376
 
377 377
   - Functions can also be used for hoppings. In this case, they take
378
-    two `~kwant.builder.Site`'s as arguments and then an arbitrary number
378
+    two `~kwant.system.Site`'s as arguments and then an arbitrary number
379 379
     of additional arguments.
380 380
 
381
-  - Apart from the real-space position `pos`, `~kwant.builder.Site` has also an
381
+  - Apart from the real-space position `pos`, `~kwant.system.Site` has also an
382 382
     attribute `tag` containing the lattice indices of the site.
383 383
 
384 384
 .. _tutorial-abring:
... ...
@@ -449,7 +449,7 @@ that returns ``True`` whenever a point is inside the shape, and
449 449
         return (r1 ** 2 < rsq < r2 ** 2)
450 450
 
451 451
 Note that this function takes a real-space position as argument (not a
452
-`~kwant.builder.Site`).
452
+`~kwant.system.Site`).
453 453
 
454 454
 We can now simply add all of the lattice points inside this shape at
455 455
 once, using the function `~kwant.lattice.Square.shape`
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2013 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -12,12 +12,16 @@ import numpy as np
12 12
 from scipy import sparse as sp
13 13
 from itertools import chain
14 14
 import types
15
+import bisect
15 16
 
16 17
 from .graph.core cimport CGraph, gintArraySlice
17 18
 from .graph.defs cimport gint
18 19
 from .graph.defs import gint_dtype
19 20
 from ._common import deprecate_args
20 21
 
22
+
23
+### Non-vectorized methods
24
+
21 25
 msg = ('Hopping from site {0} to site {1} does not match the '
22 26
        'dimensions of onsite Hamiltonians of these sites.')
23 27
 
... ...
@@ -352,3 +356,383 @@ def hamiltonian_submatrix(self, args=(), to_sites=None, from_sites=None,
352 356
         mat = func(ham, args, params, self.graph, diag, from_sites,
353 357
                    n_by_to_site, to_norb, to_off, from_norb, from_off)
354 358
     return (mat, to_norb, from_norb) if return_norb else mat
359
+
360
+
361
+### Vectorized methods
362
+
363
+
364
+_shape_error_msg = (
365
+    "The following hoppings have matrix elements of incompatible shape "
366
+    "with other matrix elements: {}"
367
+)
368
+
369
+
370
+@cython.boundscheck(False)
371
+def _vectorized_make_sparse(subgraphs, hams, long [:] norbs, long [:] orb_offsets,
372
+                 long [:] site_offsets, shape=None):
373
+    ndata = sum(h.shape[0] * h.shape[1] * h.shape[2] for h in hams)
374
+
375
+    cdef long [:] rows_view, cols_view
376
+    cdef complex [:] data_view
377
+
378
+    rows = np.empty((ndata,), dtype=long)
379
+    cols = np.empty((ndata,), dtype=long)
380
+    data = np.empty((ndata,), dtype=complex)
381
+    rows_view = rows
382
+    cols_view = cols
383
+    data_view = data
384
+
385
+    cdef long i, j, k, m, N, M, P, to_off, from_off,\
386
+              ta, fa, to_norbs, from_norbs
387
+    cdef long [:] ts_offs, fs_offs
388
+    cdef complex [:, :, :] h
389
+
390
+    m = 0
391
+    # This outer loop zip() is pure Python, but that's ok, as it
392
+    # has very few entries and the inner loops are fully vectorized
393
+    for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams):
394
+        N = h.shape[0]
395
+        M = h.shape[1]
396
+        P = h.shape[2]
397
+
398
+        if norbs[ta] != M or norbs[fa] != P:
399
+            to_sites = site_offsets[ta] + np.array(ts_offs)
400
+            from_sites = site_offsets[fa] + np.array(fs_offs)
401
+            hops = np.array([to_sites, from_sites]).transpose()
402
+            raise ValueError(_shape_error_msg.format(hops))
403
+
404
+        for i in range(N):
405
+            to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i]
406
+            from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i]
407
+            for j in range(M):
408
+                for k in range(P):
409
+                    rows_view[m] = to_off + j
410
+                    cols_view[m] = from_off + k
411
+                    data_view[m] = h[i, j, k]
412
+                    m += 1
413
+
414
+    if shape is None:
415
+        shape = (orb_offsets[-1], orb_offsets[-1])
416
+
417
+    return sp.coo_matrix((data, (rows, cols)), shape=shape)
418
+
419
+
420
+@cython.boundscheck(False)
421
+def _vectorized_make_dense(subgraphs, hams, long [:] norbs, long [:] orb_offsets,
422
+                long [:] site_offsets, shape=None):
423
+    if shape is None:
424
+        shape = (orb_offsets[-1], orb_offsets[-1])
425
+    mat = np.zeros(shape, dtype=complex)
426
+    cdef complex [:, :] mat_view
427
+    mat_view = mat
428
+
429
+    cdef long i, j, k, N, M, P, to_off, from_off,\
430
+              ta, fa, to_norbs, from_norbs
431
+    cdef long [:] ts_offs, fs_offs
432
+    cdef complex [:, :, :] h
433
+
434
+    # This outer loop zip() is pure Python, but that's ok, as it
435
+    # has very few entries and the inner loops are fully vectorized
436
+    for ((ta, fa), (ts_offs, fs_offs)), h in zip(subgraphs, hams):
437
+        N = h.shape[0]
438
+        M = h.shape[1]
439
+        P = h.shape[2]
440
+
441
+        if norbs[ta] != M or norbs[fa] != P:
442
+            to_sites = site_offsets[ta] + np.array(ts_offs)
443
+            from_sites = site_offsets[fa] + np.array(fs_offs)
444
+            hops = np.array([to_sites, from_sites]).transpose()
445
+            raise ValueError(_shape_error_msg.format(hops))
446
+
447
+        for i in range(N):
448
+            to_off = orb_offsets[ta] + norbs[ta] * ts_offs[i]
449
+            from_off = orb_offsets[fa] + norbs[fa] * fs_offs[i]
450
+            for j in range(M):
451
+                for k in range(P):
452
+                    mat_view[to_off + j, from_off + k] = h[i, j, k]
453
+    return mat
454
+
455
+
456
+def _count_norbs(syst, site_offsets, hams, args=(), params=None):
457
+    """Return the norbs and orbital offset of each site family in 'syst'
458
+
459
+    Parameters
460
+    ----------
461
+    syst : `kwant.system.System`
462
+    site_offsets : sequence of int
463
+        Contains the index of the first site for each site array
464
+        in 'syst.site_arrays'.
465
+    hams : sequence of arrays or 'None'
466
+        The Hamiltonian for each term in 'syst.terms'. 'None'
467
+        if the corresponding term has not been evaluated.
468
+    args, params : positional and keyword arguments to the system Hamiltonian
469
+    """
470
+    site_ranges = syst.site_ranges
471
+    if site_ranges is None:
472
+        # NOTE: Remove this branch once we can rely on
473
+        #       site families storing the norb information.
474
+
475
+        site_arrays = syst.site_arrays
476
+        family_norbs = {s.family: None for s in site_arrays}
477
+
478
+        # Compute the norbs per site family using already evaluated
479
+        # Hamiltonian terms where possible
480
+        for ham, t in zip(hams, syst.terms):
481
+            (to_w, from_w), _ = syst.subgraphs[t.subgraph]
482
+            if ham is not None:
483
+                family_norbs[site_arrays[to_w].family] = ham.shape[1]
484
+                family_norbs[site_arrays[from_w].family] = ham.shape[2]
485
+
486
+        # Evaluate Hamiltonian terms where we do not already have them
487
+        for n, t in enumerate(syst.terms):
488
+            (to_w, from_w), _ = syst.subgraphs[t.subgraph]
489
+            to_fam = site_arrays[to_w].family
490
+            from_fam = site_arrays[from_w].family
491
+            if family_norbs[to_fam] is None or family_norbs[from_fam] is None:
492
+                ham = syst.hamiltonian_term(n, args=args, params=params)
493
+                family_norbs[to_fam] = ham.shape[1]
494
+                family_norbs[from_fam] = ham.shape[2]
495
+
496
+        # This case is technically allowed by the format (some sites present
497
+        # but no hamiltonian terms that touch them) but is very unlikely.
498
+        if any(norbs is None for norbs in family_norbs.values()):
499
+            raise ValueError("Cannot determine the number of orbitals for "
500
+                             "some site families.")
501
+
502
+        orb_offset = 0
503
+        site_ranges = []
504
+        for start, site_array in zip(site_offsets, syst.site_arrays):
505
+            norbs = family_norbs[site_array.family]
506
+            site_ranges.append((start, norbs, orb_offset))
507
+            orb_offset += len(site_array) * norbs
508
+        site_ranges.append((site_offsets[-1], 0, orb_offset))
509
+        site_ranges = np.array(site_ranges)
510
+
511
+    _, norbs, orb_offsets = site_ranges.transpose()
512
+    # The final (extra) element in orb_offsets is the total number of orbitals
513
+    assert len(orb_offsets) == len(syst.site_arrays) + 1
514
+
515
+    return norbs, orb_offsets
516
+
517
+
518
+def _expand_norbs(compressed_norbs, site_offsets):
519
+    "Return norbs per site, given norbs per site array."
520
+    norbs = np.empty(site_offsets[-1], int)
521
+    for start, stop, norb in zip(site_offsets, site_offsets[1:],
522
+                                  compressed_norbs):
523
+        norbs[start:stop] = norb
524
+    return norbs
525
+
526
+
527
+def _reverse_subgraph(subgraph):
528
+    (to_sa, from_sa), (to_off, from_off) = subgraph
529
+    return ((from_sa, to_sa), (from_off, to_off))
530
+
531
+
532
+@deprecate_args
533
+@cython.binding(True)
534
+def vectorized_hamiltonian_submatrix(self, args=(), sparse=False,
535
+                                     return_norb=False, *, params=None):
536
+    """Return The system Hamiltonian.
537
+
538
+    Parameters
539
+    ----------
540
+    args : tuple, defaults to empty
541
+        Positional arguments to pass to ``hamiltonian_term``. Mutually
542
+        exclusive with 'params'.
543
+    sparse : bool
544
+        Whether to return a sparse or a dense matrix. Defaults to ``False``.
545
+    return_norb : bool
546
+        Whether to return arrays of numbers of orbitals.  Defaults to ``False``.
547
+    params : dict, optional
548
+        Dictionary of parameter names and their values. Mutually exclusive
549
+        with 'args'.
550
+
551
+    Returns
552
+    -------
553
+    hamiltonian_part : numpy.ndarray or scipy.sparse.coo_matrix
554
+       The Hamiltonian of the system.
555
+    norb : array of integers
556
+        Numbers of orbitals on each site. Only returned when ``return_norb``
557
+        is true.
558
+
559
+    Notes
560
+    -----
561
+    Providing positional arguments via 'args' is deprecated,
562
+    instead, provide named parameters as a dictionary via 'params'.
563
+    """
564
+
565
+    site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
566
+
567
+    subgraphs = [
568
+        self.subgraphs[t.subgraph]
569
+        for t in self.terms
570
+    ]
571
+    # Add Hermitian conjugates
572
+    subgraphs += [
573
+        _reverse_subgraph(self.subgraphs[t.subgraph])
574
+        for t in self.terms
575
+        if t.hermitian
576
+    ]
577
+
578
+    hams = [
579
+        self.hamiltonian_term(n, args=args, params=params)
580
+        for n, _ in enumerate(self.terms)
581
+    ]
582
+    # Add Hermitian conjugates
583
+    hams += [
584
+        ham.conjugate().transpose(0, 2, 1)
585
+        for ham, t in zip(hams, self.terms)
586
+        if t.hermitian
587
+    ]
588
+
589
+    norbs, orb_offsets = _count_norbs(self, site_offsets, hams,
590
+                                      args=args, params=params)
591
+
592
+    func = _vectorized_make_sparse if sparse else _vectorized_make_dense
593
+    mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets)
594
+
595
+    if return_norb:
596
+        return (mat, _expand_norbs(norbs, site_offsets))
597
+    else:
598
+        return mat
599
+
600
+
601
+@deprecate_args
602
+@cython.binding(True)
603
+def vectorized_cell_hamiltonian(self, args=(), sparse=False, *, params=None):
604
+    """Hamiltonian of a single cell of the infinite system.
605
+
606
+    Providing positional arguments via 'args' is deprecated,
607
+    instead, provide named parameters as a dictionary via 'params'.
608
+    """
609
+
610
+    site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
611
+    # Site array where next cell starts
612
+    next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
613
+
614
+    def inside_cell(term):
615
+        (to_which, from_which), _= self.subgraphs[term.subgraph]
616
+        return to_which < next_cell and from_which < next_cell
617
+
618
+    cell_terms = [
619
+        n for n, t in enumerate(self.terms)
620
+        if inside_cell(t)
621
+    ]
622
+
623
+    subgraphs = [
624
+        self.subgraphs[self.terms[n].subgraph]
625
+        for n in cell_terms
626
+    ]
627
+    # Add Hermitian conjugates
628
+    subgraphs += [
629
+        _reverse_subgraph(self.subgraphs[self.terms[n].subgraph])
630
+        for n in cell_terms
631
+        if self.terms[n].hermitian
632
+    ]
633
+
634
+    hams = [
635
+        self.hamiltonian_term(n, args=args, params=params)
636
+        for n in cell_terms
637
+    ]
638
+    # Add Hermitian conjugates
639
+    hams += [
640
+        ham.conjugate().transpose(0, 2, 1)
641
+        for ham, n in zip(hams, cell_terms)
642
+        if self.terms[n].hermitian
643
+    ]
644
+
645
+    # _count_norbs requires passing hamiltonians for all terms, or 'None'
646
+    # if it has not been evaluated
647
+    all_hams = [None] * len(self.terms)
648
+    for n, ham in zip(cell_terms, hams):
649
+        all_hams[n] = ham
650
+    norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams,
651
+                                      args=args, params=params)
652
+
653
+    shape = (orb_offsets[next_cell], orb_offsets[next_cell])
654
+    func = _vectorized_make_sparse if sparse else _vectorized_make_dense
655
+    mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
656
+
657
+    return mat
658
+
659
+
660
+@deprecate_args
661
+@cython.binding(True)
662
+def vectorized_inter_cell_hopping(self, args=(), sparse=False, *, params=None):
663
+    """Hopping Hamiltonian between two cells of the infinite system.
664
+
665
+    Providing positional arguments via 'args' is deprecated,
666
+    instead, provide named parameters as a dictionary via 'params'.
667
+    """
668
+
669
+    # Take advantage of the fact that there are distinct entries in
670
+    # onsite_terms for sites inside the unit cell, and distinct entries
671
+    # in onsite_terms for hoppings between the unit cell sites and
672
+    # interface sites. This way we only need to check the first entry
673
+    # in onsite/hopping_terms
674
+
675
+    site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
676
+    # Site array where next cell starts
677
+    next_cell = bisect.bisect(site_offsets, self.cell_size) - 1
678
+
679
+    inter_cell_hopping_terms = [
680
+        n for n, t in enumerate(self.terms)
681
+        # *from* site is in interface
682
+        if (self.subgraphs[t.subgraph][0][1] >= next_cell
683
+            and self.subgraphs[t.subgraph][0][0] < next_cell)
684
+    ]
685
+    reversed_inter_cell_hopping_terms = [
686
+        n for n, t in enumerate(self.terms)
687
+        # *to* site is in interface
688
+        if (self.subgraphs[t.subgraph][0][0] >= next_cell
689
+            and self.subgraphs[t.subgraph][0][1] < next_cell)
690
+    ]
691
+
692
+    # Evaluate inter-cell hoppings only
693
+    inter_cell_hams = [
694
+        self.hamiltonian_term(n, args=args, params=params)
695
+        for n in inter_cell_hopping_terms
696
+    ]
697
+    reversed_inter_cell_hams = [
698
+        self.hamiltonian_term(n, args=args, params=params)
699
+            .conjugate().transpose(0, 2, 1)
700
+        for n in reversed_inter_cell_hopping_terms
701
+    ]
702
+
703
+    hams = inter_cell_hams + reversed_inter_cell_hams
704
+
705
+    subgraphs = [
706
+        self.subgraphs[self.terms[n].subgraph]
707
+        for n in inter_cell_hopping_terms
708
+    ]
709
+    subgraphs += [
710
+        _reverse_subgraph(self.subgraphs[self.terms[n].subgraph])
711
+        for n in reversed_inter_cell_hopping_terms
712
+    ]
713
+
714
+    # All the 'from' sites are in the previous domain, but to build a
715
+    # matrix we need to get the equivalent sites in the fundamental domain.
716
+    # We set the 'from' site array to the one from the fundamental domain.
717
+    subgraphs = [
718
+        ((to_sa, from_sa - next_cell), (to_off, from_off))
719
+        for (to_sa, from_sa), (to_off, from_off) in subgraphs
720
+    ]
721
+
722
+    # _count_norbs requires passing hamiltonians for all terms, or 'None'
723
+    # if it has not been evaluated
724
+    all_hams = [None] * len(self.terms)
725
+    for n, ham in zip(inter_cell_hopping_terms, inter_cell_hams):
726
+        all_hams[n] = ham
727
+    for n, ham in zip(reversed_inter_cell_hopping_terms,
728
+                      reversed_inter_cell_hams):
729
+        # Transpose to get back correct shape wrt. original term
730
+        all_hams[n] = ham.transpose(0, 2, 1)
731
+    norbs, orb_offsets = _count_norbs(self, site_offsets, all_hams,
732
+                                      args=args, params=params)
733
+
734
+    shape = (orb_offsets[next_cell],
735
+             orb_offsets[len(self.site_arrays) - next_cell])
736
+    func = _vectorized_make_sparse if sparse else _vectorized_make_dense
737
+    mat = func(subgraphs, hams, norbs, orb_offsets, site_offsets, shape=shape)
738
+    return mat
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2016 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -14,11 +14,14 @@ import copy
14 14
 from functools import total_ordering, wraps, update_wrapper
15 15
 from itertools import islice, chain
16 16
 import textwrap
17
+import bisect
18
+import numbers
17 19
 import inspect
18 20
 import tinyarray as ta
19 21
 import numpy as np
20 22
 from scipy import sparse
21 23
 from . import system, graph, KwantDeprecationWarning, UserCodeError
24
+from .system import Site, SiteArray, SiteFamily
22 25
 from .linalg import lll
23 26
 from .operator import Density
24 27
 from .physics import DiscreteSymmetry, magnetic_gauge
... ...
@@ -26,182 +29,13 @@ from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
26 29
                       interleave, deprecate_args, memoize)
27 30
 
28 31
 
29
-__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
30
-           'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead',
31
-           'add_peierls_phase']
32
+__all__ = ['Builder', 'SimpleSiteFamily', 'Symmetry', 'HoppingKind', 'Lead',
33
+           'BuilderLead', 'SelfEnergyLead', 'ModesLead', 'add_peierls_phase']
32 34
 
33 35
 
34
-################ Sites and site families
35
-
36
-class Site(tuple):
37
-    """A site, member of a `SiteFamily`.
38
-
39
-    Sites are the vertices of the graph which describes the tight binding
40
-    system in a `Builder`.
41
-
42
-    A site is uniquely identified by its family and its tag.
43
-
44
-    Parameters
45
-    ----------
46
-    family : an instance of `SiteFamily`
47
-        The 'type' of the site.
48
-    tag : a hashable python object
49
-        The unique identifier of the site within the site family, typically a
50
-        vector of integers.
51
-
52
-    Raises
53
-    ------
54
-    ValueError
55
-        If `tag` is not a proper tag for `family`.
56
-
57
-    Notes
58
-    -----
59
-    For convenience, ``family(*tag)`` can be used instead of ``Site(family,
60
-    tag)`` to create a site.
61
-
62
-    The parameters of the constructor (see above) are stored as instance
63
-    variables under the same names.  Given a site ``site``, common things to
64
-    query are thus ``site.family``, ``site.tag``, and ``site.pos``.
65
-    """
66
-    __slots__ = ()
67
-
68
-    family = property(operator.itemgetter(0),
69
-                      doc="The site family to which the site belongs.")
70
-    tag = property(operator.itemgetter(1), doc="The tag of the site.")
71
-
72
-
73
-    def __new__(cls, family, tag, _i_know_what_i_do=False):
74
-        if _i_know_what_i_do:
75
-            return tuple.__new__(cls, (family, tag))
76
-        try:
77
-            tag = family.normalize_tag(tag)
78
-        except (TypeError, ValueError) as e:
79
-            msg = 'Tag {0} is not allowed for site family {1}: {2}'
80
-            raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
81
-        return tuple.__new__(cls, (family, tag))
82
-
83
-    def __repr__(self):
84
-        return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
85
-
86
-    def __str__(self):
87
-        sf = self.family
88
-        return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
89
-
90
-    def __getnewargs__(self):
91
-        return (self.family, self.tag, True)
92
-
93
-    @property
94
-    def pos(self):
95
-        """Real space position of the site.
96
-
97
-        This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
98
-        """
99
-        return self.family.pos(self.tag)
100
-
36
+################ Site families
101 37
 
102 38
 @total_ordering
103
-class SiteFamily(metaclass=abc.ABCMeta):
104
-    """Abstract base class for site families.
105
-
106
-    Site families are the 'type' of `Site` objects.  Within a family, individual
107
-    sites are uniquely identified by tags.  Valid tags must be hashable Python
108
-    objects, further details are up to the family.
109
-
110
-    Site families must be immutable and fully defined by their initial
111
-    arguments.  They must inherit from this abstract base class and call its
112
-    __init__ function providing it with two arguments: a canonical
113
-    representation and a name.  The canonical representation will be returned as
114
-    the objects representation and must uniquely identify the site family
115
-    instance.  The name is a string used to distinguish otherwise identical site
116
-    families.  It may be empty. ``norbs`` defines the number of orbitals
117
-    on sites associated with this site family; it may be `None`, in which case
118
-    the number of orbitals is not specified.
119
-
120
-
121
-    All site families must define the method `normalize_tag` which brings a tag
122
-    to the standard format for this site family.
123
-
124
-    Site families that are intended for use with plotting should also provide a
125
-    method `pos(tag)`, which returns a vector with real-space coordinates of the
126
-    site belonging to this family with a given tag.
127
-
128
-    If the ``norbs`` of a site family are provided, and sites of this family
129
-    are used to populate a `~kwant.builder.Builder`, then the associated
130
-    Hamiltonian values must have the correct shape. That is, if a site family
131
-    has ``norbs = 2``, then any on-site terms for sites belonging to this
132
-    family should be 2x2 matrices. Similarly, any hoppings to/from sites
133
-    belonging to this family must have a matrix structure where there are two
134
-    rows/columns. This condition applies equally to Hamiltonian values that
135
-    are given by functions. If this condition is not satisfied, an error will
136
-    be raised.
137
-    """
138
-
139
-    def __init__(self, canonical_repr, name, norbs):
140
-        self.canonical_repr = canonical_repr
141
-        self.hash = hash(canonical_repr)
142
-        self.name = name
143
-        if norbs is None:
144
-            warnings.warn("Not specfying norbs is deprecated. Always specify "
145
-                          "norbs when creating site families.",
146
-                          KwantDeprecationWarning, stacklevel=3)
147
-        if norbs is not None:
148
-            if int(norbs) != norbs or norbs <= 0:
149
-                raise ValueError('The norbs parameter must be an integer > 0.')
150
-            norbs = int(norbs)
151
-        self.norbs = norbs
152
-
153
-    def __repr__(self):
154
-        return self.canonical_repr
155
-
156
-    def __str__(self):
157
-        if self.name:
158
-            msg = '<{0} site family {1}{2}>'
159
-        else:
160
-            msg = '<unnamed {0} site family{2}>'
161
-        orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
162
-        return msg.format(self.__class__.__name__, self.name, orbs)
163
-
164
-    def __hash__(self):
165
-        return self.hash
166
-
167
-    def __eq__(self, other):
168
-        try:
169
-            return self.canonical_repr == other.canonical_repr
170
-        except AttributeError:
171
-            return False
172
-
173
-    def __ne__(self, other):
174
-        try:
175
-            return self.canonical_repr != other.canonical_repr
176
-        except AttributeError:
177
-            return True
178
-
179
-    def __lt__(self, other):
180
-        # If this raises an AttributeError, we were trying
181
-        # to compare it to something non-comparable anyway.
182
-        return self.canonical_repr < other.canonical_repr
183
-
184
-    @abc.abstractmethod
185
-    def normalize_tag(self, tag):
186
-        """Return a normalized version of the tag.
187
-
188
-        Raises TypeError or ValueError if the tag is not acceptable.
189
-        """
190
-        pass
191
-
192
-    def __call__(self, *tag):
193
-        """
194
-        A convenience function.
195
-
196
-        This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
197
-        """
198
-        # Catch a likely and difficult to find mistake.
199
-        if tag and isinstance(tag[0], tuple):
200
-            raise ValueError('Use site_family(1, 2) instead of '
201
-                             'site_family((1, 2))!')
202
-        return Site(self, tag)
203
-
204
-
205 39
 class SimpleSiteFamily(SiteFamily):
206 40
     """A site family used as an example and for testing.
207 41
 
... ...
@@ -307,19 +141,44 @@ class Symmetry(metaclass=abc.ABCMeta):
307 141
     def which(self, site):
308 142
         """Calculate the domain of the site.
309 143
 
310
-        Return the group element whose action on a certain site from the
311
-        fundamental domain will result in the given ``site``.
144
+        Parameters
145
+        ----------
146
+        site : `~kwant.system.Site` or `~kwant.system.SiteArray`
147
+
148
+        Returns
149
+        -------
150
+        group_element : tuple or sequence of tuples
151
+            A single tuple if ``site`` is a Site, or a sequence of tuples if
152
+            ``site`` is a SiteArray.  The group element(s) whose action
153
+            on a certain site(s) from the fundamental domain will result
154
+            in the given ``site``.
312 155
         """
313 156
         pass
314 157
 
315 158
     @abc.abstractmethod
316 159
     def act(self, element, a, b=None):
317
-        """Act with a symmetry group element on a site or hopping."""
160
+        """Act with symmetry group element(s) on site(s) or hopping(s).
161
+
162
+        Parameters
163
+        ----------
164
+        element : tuple or sequence of tuples
165
+            Group element(s) with which to act on the provided site(s)
166
+            or hopping(s)
167
+        a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
168
+            If Site then ``element`` is a single tuple, if SiteArray then
169
+            ``element`` is a sequence of tuples. If only ``a`` is provided then
170
+            ``element`` acts on the site(s) of ``a``. If ``b`` is also provided
171
+            then ``element`` acts on the hopping(s) ``(a, b)``.
172
+        """
318 173
         pass
319 174
 
320 175
     def to_fd(self, a, b=None):
321 176
         """Map a site or hopping to the fundamental domain.
322 177
 
178
+        Parameters
179
+        ----------
180
+        a, b : `~kwant.system.Site` or `~kwant.system.SiteArray`
181
+
323 182
         If ``b`` is None, return a site equivalent to ``a`` within the
324 183
         fundamental domain.  Otherwise, return a hopping equivalent to ``(a,
325 184
         b)`` but where the first element belongs to the fundamental domain.
... ...
@@ -329,11 +188,30 @@ class Symmetry(metaclass=abc.ABCMeta):
329 188
         return self.act(-self.which(a), a, b)
330 189
 
331 190
     def in_fd(self, site):
332
-        """Tell whether ``site`` lies within the fundamental domain."""
333
-        for d in self.which(site):
334
-            if d != 0:
335
-                return False
336
-        return True
191
+        """Tell whether ``site`` lies within the fundamental domain.
192
+
193
+        Parameters
194
+        ----------
195
+        site : `~kwant.system.Site` or `~kwant.system.SiteArray`
196
+
197
+        Returns
198
+        -------
199
+        in_fd : bool or sequence of bool
200
+            single bool if ``site`` is a Site, or a sequence of
201
+            bool if ``site`` is a SiteArray. In the latter case
202
+            we return whether each site in the SiteArray is in
203
+            the fundamental domain.
204
+        """
205
+        if isinstance(site, Site):
206
+            for d in self.which(site):
207
+                if d != 0:
208
+                    return False
209
+            return True
210
+        elif isinstance(site, SiteArray):
211
+            which = self.which(site)
212
+            return np.logical_and.reduce(which != 0, axis=1)
213
+        else:
214
+            raise TypeError("'site' must be a Site or SiteArray")
337 215
 
338 216
     @abc.abstractmethod
339 217
     def subgroup(self, *generators):
... ...
@@ -412,8 +290,8 @@ class HoppingKind(tuple):
412 290
     ----------
413 291
     delta : Sequence of integers
414 292
         The sequence is interpreted as a vector with integer elements.
415
-    family_a : `~kwant.builder.SiteFamily`
416
-    family_b : `~kwant.builder.SiteFamily` or ``None`` (default)
293
+    family_a : `~kwant.system.SiteFamily`
294
+    family_b : `~kwant.system.SiteFamily` or ``None`` (default)
417 295
         The default value means: use the same family as `family_a`.
418 296
 
419 297
     Notes
... ...
@@ -570,7 +448,7 @@ class BuilderLead(Lead):
570 448
         The tight-binding system of a lead. It has to possess appropriate
571 449
         symmetry, and it may not contain hoppings between further than
572 450
         neighboring images of the fundamental domain.
573
-    interface : sequence of `Site` instances
451
+    interface : sequence of `~kwant.system.Site` instances
574 452
         Sequence of sites in the scattering region to which the lead is
575 453
         attached.
576 454
 
... ...
@@ -578,9 +456,9 @@ class BuilderLead(Lead):
578 456
     ----------
579 457
     builder : `Builder`
580 458
         The tight-binding system of a lead.
581
-    interface : list of `Site` instances
459
+    interface : list of `~kwant.system.Site` instances
582 460
         A sorted list of interface sites.
583
-    padding : list of `Site` instances
461
+    padding : list of `~kwant.system.Site` instances
584 462
         A sorted list of sites that originate from the lead, have the same
585 463
         onsite Hamiltonian, and are connected by the same hoppings as the lead
586 464
         sites.
... ...
@@ -607,7 +485,10 @@ class BuilderLead(Lead):
607 485
 
608 486
         The order of interface sites is kept during finalization.
609 487
         """
610
-        return InfiniteSystem(self.builder, self.interface)
488
+        if self.builder.vectorize:
489
+            return InfiniteVectorizedSystem(self.builder, self.interface)
490
+        else:
491
+            return InfiniteSystem(self.builder, self.interface)
611 492
 
612 493
 
613 494
 def _ensure_signature(func):
... ...
@@ -638,7 +519,7 @@ class SelfEnergyLead(Lead):
638 519
     selfenergy_func : function
639 520
         Has the same signature as `selfenergy` (without the ``self``
640 521
         parameter) and returns the self energy matrix for the interface sites.
641
-    interface : sequence of `Site` instances
522
+    interface : sequence of `~kwant.system.Site` instances
642 523
     parameters : sequence of strings
643 524
         The parameters on which the lead depends.
644 525
     """
... ...
@@ -669,7 +550,7 @@ class ModesLead(Lead):
669 550
         and returns the modes of the lead as a tuple of
670 551
         `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`.
671 552
     interface :
672
-        sequence of `Site` instances
553
+        sequence of `~kwant.system.Site` instances
673 554
     parameters : sequence of strings
674 555
         The parameters on which the lead depends.
675 556
     """
... ...
@@ -797,10 +678,11 @@ class Builder:
797 678
     This is one of the central types in Kwant.  It is used to construct tight
798 679
     binding systems in a flexible way.
799 680
 
800
-    The nodes of the graph are `Site` instances.  The edges, i.e. the hoppings,
801
-    are pairs (2-tuples) of sites.  Each node and each edge has a value
802
-    associated with it.  The values associated with nodes are interpreted as
803
-    on-site Hamiltonians, the ones associated with edges as hopping integrals.
681
+    The nodes of the graph are `~kwant.system.Site` instances.  The edges,
682
+    i.e. the hoppings, are pairs (2-tuples) of sites.  Each node and each
683
+    edge has a value associated with it.  The values associated with nodes
684
+    are interpreted as on-site Hamiltonians, the ones associated with edges
685
+    as hopping integrals.
804 686
 
805 687
     To make the graph accessible in a way that is natural within the Python
806 688
     language it is exposed as a *mapping* (much like a built-in Python
... ...
@@ -827,6 +709,14 @@ class Builder:
827 709
     chiral : 2D array, dictionary, function or `None`
828 710
         The unitary part of the onsite chiral symmetry operator.
829 711
         Same format as that of `conservation_law`.
712
+    vectorize : bool, default: False
713
+        If True then the finalized Builder will evaluate its Hamiltonian in
714
+        a vectorized manner. This requires that all value functions take
715
+        `~kwant.system.SiteArray` as first/second parameter rather than
716
+        `~kwant.system.Site`, and return an array of values. The returned
717
+        array must have the same length as the provided SiteArray(s), and
718
+        may contain either scalars (i.e. a vector of values) or matrices
719
+        (i.e. a 3d array of values).
830 720
 
831 721
     Notes
832 722
     -----
... ...
@@ -905,7 +795,7 @@ class Builder:
905 795
     """
906 796
 
907 797
     def __init__(self, symmetry=None, *, conservation_law=None, time_reversal=None,
908
-                 particle_hole=None, chiral=None):
798
+                 particle_hole=None, chiral=None, vectorize=False):
909 799
         if symmetry is None:
910 800
             symmetry = NoSymmetry()
911 801
         else:
... ...
@@ -915,6 +805,7 @@ class Builder:
915 805
         self.time_reversal = time_reversal
916 806
         self.particle_hole = particle_hole
917 807
         self.chiral = chiral
808
+        self.vectorize = vectorize
918 809
         self.leads = []
919 810
         self.H = {}
920 811
 
... ...
@@ -998,6 +889,7 @@ class Builder:
998 889
         result.time_reversal = self.time_reversal
999 890
         result.particle_hole = self.particle_hole
1000 891
         result.chiral = self.chiral
892
+        result.vectorize = self.vectorize
1001 893
         result.leads = self.leads
1002 894
         result.H = self.H
1003 895
         return result
... ...
@@ -1452,7 +1344,7 @@ class Builder:
1452 1344
             A boolean function of site returning whether the site should be
1453 1345
             added to the target builder or not. The shape must be compatible
1454 1346
             with the symmetry of the target builder.
1455
-        start : `Site` instance or iterable thereof or iterable of numbers
1347
+        start : `~kwant.system.Site` or iterable thereof or iterable of numbers
1456 1348
             The site(s) at which the the flood-fill starts.  If start is an
1457 1349
             iterable of numbers, the starting site will be
1458 1350
             ``template.closest(start)``.
... ...
@@ -1462,7 +1354,7 @@ class Builder:
1462 1354
 
1463 1355
         Returns
1464 1356
         -------
1465
-        added_sites : list of `Site` objects that were added to the system.
1357
+        added_sites : list of `~kwant.system.Site` that were added to the system.
1466 1358
 
1467 1359
         Raises
1468 1360
         ------
... ...
@@ -1614,7 +1506,7 @@ class Builder:
1614 1506
         ----------
1615 1507
         lead_builder : `Builder` with 1D translational symmetry
1616 1508
             Builder of the lead which has to be attached.
1617
-        origin : `Site`
1509
+        origin : `~kwant.system.Site`
1618 1510
             The site which should belong to a domain where the lead should
1619 1511
             begin. It is used to attach a lead inside the system, e.g. to an
1620 1512
             inner radius of a ring.
... ...
@@ -1624,7 +1516,7 @@ class Builder:
1624 1516
 
1625 1517
         Returns
1626 1518
         -------
1627
-        added_sites : list of `Site` objects that were added to the system.
1519
+        added_sites : list of `~kwant.system.Site`
1628 1520
 
1629 1521
         Raises
1630 1522
         ------
... ...
@@ -1677,6 +1569,13 @@ class Builder:
1677 1569
         if add_cells < 0 or int(add_cells) != add_cells:
1678 1570
             raise ValueError('add_cells must be an integer >= 0.')
1679 1571
 
1572
+        if self.vectorize != lead_builder.vectorize:
1573
+            raise ValueError(
1574
+                "Sites of the lead were added to the scattering "
1575
+                "region, but only one of these systems is vectorized. "
1576
+                "Vectorize both systems to allow attaching leads."
1577
+            )
1578
+
1680 1579
         sym = lead_builder.symmetry
1681 1580
         H = lead_builder.H
1682 1581
 
... ...
@@ -1697,6 +1596,7 @@ class Builder:
1697 1596
                 time_reversal=lead_builder.time_reversal,
1698 1597
                 particle_hole=lead_builder.particle_hole,
1699 1598
                 chiral=lead_builder.chiral,
1599
+                vectorize=lead_builder.vectorize,
1700 1600
             )
1701 1601
             with reraise_warnings():
1702 1602
                 new_lead.fill(lead_builder, lambda site: True,
... ...
@@ -1793,9 +1693,15 @@ class Builder:
1793 1693
         `Symmetry` can be finalized.
1794 1694
         """
1795 1695
         if self.symmetry.num_directions == 0:
1796
-            return FiniteSystem(self)
1696
+            if self.vectorize:
1697
+                return FiniteVectorizedSystem(self)
1698
+            else:
1699
+                return FiniteSystem(self)
1797 1700
         elif self.symmetry.num_directions == 1:
1798
-            return InfiniteSystem(self)
1701
+            if self.vectorize:
1702
+                return InfiniteVectorizedSystem(self)
1703
+            else:
1704
+                return InfiniteSystem(self)
1799 1705
         else:
1800 1706
             raise ValueError('Currently, only builders without or with a 1D '
1801 1707
                              'translational symmetry can be finalized.')
... ...
@@ -1969,6 +1875,10 @@ def add_peierls_phase(syst, peierls_parameter='phi', fix_gauge=True):
1969 1875
 
1970 1876
         return f
1971 1877
 
1878
+    if syst.vectorize:
1879
+        raise TypeError("'add_peierls_phase' does not work with "
1880
+                        "vectorized Builders")
1881
+
1972 1882
     ret = _add_peierls_phase(syst, peierls_parameter).finalized()
1973 1883
 
1974 1884
     if fix_gauge:
... ...
@@ -2115,6 +2025,135 @@ class _FinalizedBuilderMixin:
2115 2025
         return DiscreteSymmetry(projectors, *(evaluate(symm) for symm in
2116 2026
                                               self._symmetries))
2117 2027
 
2028
+    def pos(self, i):
2029
+        return self.sites[i].pos
2030
+
2031
+
2032
+class _VectorizedFinalizedBuilderMixin(_FinalizedBuilderMixin):
2033
+    """Common functionality for all vectorized finalized builders
2034
+
2035
+    Attributes
2036
+    ----------
2037
+    _term_values : sequence
2038
+        Each item is either an array of values (if 'param_names is None')
2039
+        or a vectorized value function (in which case 'param_names' is a
2040
+        sequence of strings: the extra parameters to the value function).
2041
+    _onsite_term_by_site_id : sequence of int
2042
+        Maps site (integer) to the term that evaluates the associated
2043
+        onsite matrix element.
2044
+    _hopping_term_by_edge_id : sequence of int
2045
+        Maps edge id in the graph to the term that evaluates the associated
2046
+        hopping matrix element. If negative, then the associated hopping is
2047
+        the Hermitian conjugate of another hopping; if the number stored is
2048
+        't' (< 0) then the associated hopping is stored in term '-t - 1'.
2049
+    """
2050
+    def hamiltonian(self, i, j, *args, params=None):
2051
+        warnings.warn(
2052
+            "Querying individual matrix elements with 'hamiltonian' is "
2053
+            "deprecated, and now takes O(log N) time where N is the number "
2054
+            "of matrix elements per hamiltonian term. Query many matrix "
2055
+            "elements at once using 'hamiltonian_term'.",
2056
+            KwantDeprecationWarning
2057
+        )
2058
+        site_offsets = np.cumsum([0] + [len(s) for s in self.site_arrays])
2059
+        if i == j:
2060
+            which_term = self._onsite_term_by_site_id[i]
2061
+            (w, _), (off, _) = self.subgraphs[self.terms[which_term].subgraph]
2062
+            i_off = i - site_offsets[w]
2063
+            selector = np.searchsorted(off, i_off)
2064
+            onsite = self.hamiltonian_term(
2065
+                which_term, [selector], args, params=params)
2066
+            return onsite[0]
2067
+        else:
2068
+            edge_id = self.graph.first_edge_id(i, j)
2069
+            which_term = self._hopping_term_by_edge_id[edge_id]
2070
+            herm_conj = which_term < 0
2071
+            if herm_conj:
2072
+                which_term = -which_term - 1
2073
+            term = self.terms[which_term]
2074
+            # To search for hopping (i, j) in hopping_terms, we need
2075
+            # to treat hopping_terms as a record array of integer pairs
2076
+            dtype = np.dtype([('f0', int), ('f1', int)])
2077
+            # For hermitian conjugate terms search through the
2078
+            # *other* term's hoppings because those are sorted.
2079
+
2080
+            (to_w, from_w), hoppings = self.subgraphs[term.subgraph]
2081
+            hoppings = hoppings.transpose()  # to get array-of-pairs
2082
+            if herm_conj:
2083
+                hop = (j - site_offsets[to_w], i - site_offsets[from_w])
2084
+            else:
2085
+                hop = (i - site_offsets[to_w], j - site_offsets[from_w])
2086
+            hop = np.array(hop, dtype=dtype)
2087
+            hoppings = hoppings.view(dtype).reshape(-1)
2088
+            selector = np.recarray.searchsorted(hoppings, hop)
2089
+            h = self.hamiltonian_term(
2090
+                    which_term, [selector], args, params=params)
2091
+            h = h[0]
2092
+            if herm_conj:
2093
+                h = h.conjugate().transpose()
2094
+            return h
2095
+
2096
+    def hamiltonian_term(self, term_number, selector=slice(None),
2097
+                         args=(), params=None):
2098
+        if args and params:
2099
+            raise TypeError("'args' and 'params' are mutually exclusive.")
2100
+
2101
+        term = self.terms[term_number]
2102
+        val = self._term_values[term_number]
2103
+
2104
+        if not callable(val):
2105
+            return val[selector]
2106
+
2107
+        # Construct site arrays to pass to the vectorized value function.
2108
+        subgraph = self.subgraphs[self.terms[term_number].subgraph]
2109
+        (to_which, from_which), (to_off, from_off) = subgraph
2110
+        is_onsite = to_off is from_off
2111
+        to_off = to_off[selector]
2112
+        from_off = from_off[selector]
2113
+        assert len(to_off) == len(from_off)
2114
+
2115
+        to_family = self.site_arrays[to_which].family
2116
+        to_tags = self.site_arrays[to_which].tags
2117
+        to_site_array = SiteArray(to_family, to_tags[to_off])
2118
+        if not is_onsite:
2119
+            from_family = self.site_arrays[from_which].family
2120
+            from_tags = self.site_arrays[from_which].tags
2121
+            from_site_array = SiteArray(from_family, from_tags[from_off])
2122
+
2123
+        # Construct args from params
2124
+        if params:
2125
+            # There was a problem extracting parameter names from the value
2126
+            # function (probably an illegal signature) and we are using
2127
+            # keyword parameters.
2128
+            if self._term_errors[term_number] is not None:
2129
+                raise self._term_errors[term_number]
2130
+            try:
2131
+                args = [params[p] for p in term.parameters]
2132
+            except KeyError:
2133
+                missing = [p for p in term.parameters if p not in params]
2134
+                msg = ('System is missing required arguments: ',
2135
+                       ', '.join(map('"{}"'.format, missing)))
2136
+                raise TypeError(''.join(msg))
2137
+
2138
+        if is_onsite:
2139
+            try:
2140
+                ham = val(to_site_array, *args)
2141
+            except Exception as exc:
2142
+                _raise_user_error(exc, val)
2143
+        else:
2144
+            try:
2145
+                ham = val(
2146
+                        *self.symmetry.to_fd(
2147
+                            to_site_array,
2148
+                            from_site_array),
2149
+                        *args)
2150
+            except Exception as exc:
2151
+                _raise_user_error(exc, val)
2152
+
2153
+        ham = _normalize_term_value(ham, len(to_site_array))
2154
+
2155
+        return ham
2156
+
2118 2157
 
2119 2158
 # The same (value, parameters) pair will be used for many sites/hoppings,
2120 2159
 # so we cache it to avoid wasting extra memory.
... ...
@@ -2147,6 +2186,64 @@ def _value_params_pair_cache(nstrip):
2147 2186
     return get
2148 2187
 
2149 2188
 
2189
+def _make_graph(H, id_by_site):
2190
+    g = graph.Graph()
2191
+    g.num_nodes = len(id_by_site)  # Some sites could not appear in any edge.
2192
+    for tail, hvhv in H.items():
2193
+        for head in islice(hvhv, 2, None, 2):
2194
+            if tail == head:
2195
+                continue
2196
+            g.add_edge(id_by_site[tail], id_by_site[head])
2197
+    return g.compressed()
2198
+
2199
+
2200
+def _finalize_leads(leads, id_by_site):
2201
+    #### Connect leads.
2202
+    finalized_leads = []
2203
+    lead_interfaces = []
2204
+    lead_paddings = []
2205
+    for lead_nr, lead in enumerate(leads):
2206
+        try:
2207
+            with warnings.catch_warnings(record=True) as ws:
2208
+                warnings.simplefilter("always")
2209
+                # The following line is the whole "payload" of the entire
2210
+                # try-block.
2211
+                finalized_leads.append(lead.finalized())
2212
+        except ValueError as e:
2213
+            # Re-raise the exception with an additional message.
2214
+            msg = 'Problem finalizing lead {0}:'.format(lead_nr)
2215
+            e.args = (' '.join((msg,) + e.args),)
2216
+            raise
2217
+        else:
2218
+            for w in ws:
2219
+                # Re-raise any warnings with an additional message and the
2220
+                # proper stacklevel.
2221
+                w = w.message
2222
+                msg = 'When finalizing lead {0}:'.format(lead_nr)
2223
+                warnings.warn(w.__class__(' '.join((msg,) + w.args)),
2224
+                              stacklevel=3)
2225
+        try:
2226
+            interface = [id_by_site[isite] for isite in lead.interface]
2227
+        except KeyError as e:
2228
+            msg = ("Lead {0} is attached to a site that does not "
2229
+                   "belong to the scattering region:\n {1}")
2230
+            raise ValueError(msg.format(lead_nr, e.args[0]))
2231
+
2232
+        lead_interfaces.append(np.array(interface))
2233
+
2234
+        padding = getattr(lead, 'padding', [])
2235
+        # Some padding sites might have been removed after the lead was
2236
+        # attached. Unlike in the case of the interface, this is not a
2237
+        # problem.
2238
+        finalized_padding = [
2239
+            id_by_site[isite] for isite in padding if isite in id_by_site
2240
+        ]
2241
+
2242
+        lead_paddings.append(np.array(finalized_padding))
2243
+
2244
+    return finalized_leads, lead_interfaces, lead_paddings
2245
+
2246
+
2150 2247
 class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2151 2248
     """Finalized `Builder` with leads.
2152 2249
 
... ...
@@ -2155,11 +2252,11 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2155 2252
     Attributes
2156 2253
     ----------
2157 2254
     sites : sequence
2158
-        ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
2255
+        ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2159 2256
         to the integer-labeled site ``i`` of the low-level system. The sites
2160 2257
         are ordered first by their family and then by their tag.
2161
-    id_by_site : dict
2162
-        The inverse of ``sites``; maps high-level `~kwant.builder.Site`
2258
+    id_by_site : mapping
2259
+        The inverse of ``sites``; maps high-level `~kwant.system.Site`
2163 2260
         instances to their integer label.
2164 2261
         Satisfies ``id_by_site[sites[i]] == i``.
2165 2262
     """
... ...
@@ -2173,57 +2270,10 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2173 2270
         for site_id, site in enumerate(sites):
2174 2271
             id_by_site[site] = site_id
2175 2272
 
2176
-        #### Make graph.
2177
-        g = graph.Graph()
2178
-        g.num_nodes = len(sites)  # Some sites could not appear in any edge.
2179
-        for tail, hvhv in builder.H.items():
2180
-            for head in islice(hvhv, 2, None, 2):
2181
-                if tail == head:
2182
-                    continue
2183
-                g.add_edge(id_by_site[tail], id_by_site[head])
2184
-        g = g.compressed()
2185
-
2186
-        #### Connect leads.
2187
-        finalized_leads = []
2188
-        lead_interfaces = []
2189
-        lead_paddings = []
2190
-        for lead_nr, lead in enumerate(builder.leads):
2191
-            try:
2192
-                with warnings.catch_warnings(record=True) as ws:
2193
-                    warnings.simplefilter("always")
2194
-                    # The following line is the whole "payload" of the entire
2195
-                    # try-block.
2196
-                    finalized_leads.append(lead.finalized())
2197
-                for w in ws:
2198
-                    # Re-raise any warnings with an additional message and the
2199
-                    # proper stacklevel.
2200
-                    w = w.message
2201
-                    msg = 'When finalizing lead {0}:'.format(lead_nr)
2202
-                    warnings.warn(w.__class__(' '.join((msg,) + w.args)),
2203
-                                  stacklevel=3)
2204
-            except ValueError as e:
2205
-                # Re-raise the exception with an additional message.
2206
-                msg = 'Problem finalizing lead {0}:'.format(lead_nr)
2207
-                e.args = (' '.join((msg,) + e.args),)
2208
-                raise
2209
-            try:
2210
-                interface = [id_by_site[isite] for isite in lead.interface]
2211
-            except KeyError as e:
2212
-                msg = ("Lead {0} is attached to a site that does not "
2213
-                       "belong to the scattering region:\n {1}")
2214
-                raise ValueError(msg.format(lead_nr, e.args[0]))
2215
-
2216
-            lead_interfaces.append(np.array(interface))
2217
-
2218
-            padding = getattr(lead, 'padding', [])
2219
-            # Some padding sites might have been removed after the lead was
2220
-            # attached. Unlike in the case of the interface, this is not a
2221
-            # problem.
2222
-            finalized_padding = [
2223
-                id_by_site[isite] for isite in padding if isite in id_by_site
2224
-            ]
2273
+        graph = _make_graph(builder.H, id_by_site)
2225 2274
 
2226
-            lead_paddings.append(np.array(finalized_padding))
2275
+        finalized_leads, lead_interfaces, lead_paddings =\
2276
+            _finalize_leads(builder.leads, id_by_site)
2227 2277
 
2228 2278
         # Because many onsites/hoppings share the same (value, parameter)
2229 2279
         # pairs, we keep them in a cache so that we only store a given pair
... ...
@@ -2232,7 +2282,7 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2232 2282
         onsites = [cache(builder.H[site][1]) for site in sites]
2233 2283
         cache = _value_params_pair_cache(2)
2234 2284
         hoppings = [cache(builder._get_edge(sites[tail], sites[head]))
2235
-                    for tail, head in g]
2285
+                    for tail, head in graph]
2236 2286
 
2237 2287
         # Compute the union of the parameters of onsites and hoppings.  Here,
2238 2288
         # 'onsites' and 'hoppings' are pairs whose second element is one of
... ...
@@ -2252,7 +2302,7 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2252 2302
         else:
2253 2303
             parameters = frozenset(parameters)
2254 2304
 
2255
-        self.graph = g
2305
+        self.graph = graph
2256 2306
         self.sites = sites
2257 2307
         self.site_ranges = _site_ranges(sites)
2258 2308
         self.id_by_site = id_by_site
... ...
@@ -2265,8 +2315,502 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
2265 2315
         self.lead_paddings = lead_paddings
2266 2316
         self._init_discrete_symmetries(builder)
2267 2317
 
2268
-    def pos(self, i):
2269
-        return self.sites[i].pos
2318
+
2319
+def _make_site_arrays(sites):
2320
+    tags_by_family = {}
2321
+    for family, tag in sites:  # Sites are tuples
2322
+        tags_by_family.setdefault(family, []).append(tag)
2323
+    site_arrays = []
2324
+    for family, tags in sorted(tags_by_family.items()):
2325
+        site_arrays.append(SiteArray(family, sorted(tags)))
2326
+    return site_arrays
2327
+
2328
+
2329
+# Wrapper for accessing sites by their sequential number from a
2330
+# list of SiteArrays.
2331
+class _Sites:
2332
+
2333
+    def __init__(self, site_arrays):
2334
+        self._site_arrays = site_arrays
2335
+        lengths = [0] + [len(arr.tags) for arr in site_arrays]
2336
+        self._start_idxs = np.cumsum(lengths)[:-1]
2337
+
2338
+    def __len__(self):
2339
+        return sum(len(arr.tags) for arr in self._site_arrays)
2340
+
2341
+    def __getitem__(self, idx):
2342
+        if not isinstance(idx, numbers.Integral):
2343
+            raise TypeError("Only individual sites may be indexed")
2344
+        if idx < 0 or idx >= len(self):
2345
+            raise IndexError(idx)
2346
+
2347
+        which = bisect.bisect(self._start_idxs, idx) - 1
2348
+        arr = self._site_arrays[which]
2349
+        start = self._start_idxs[which]
2350
+        fam = arr.family
2351
+        tag = arr.tags[idx - start]
2352
+        return Site(fam, tag)
2353
+
2354
+    def __iter__(self):
2355
+        for arr in self._site_arrays:
2356
+            for tag in arr.tags:
2357
+                yield Site(arr.family, tag)
2358
+
2359
+
2360
+# Wrapper for accessing the sequential number of a site, given
2361
+# Site object, from a list of SiteArrays.
2362
+class _IdBySite:
2363
+
2364
+    def __init__(self, site_arrays):
2365
+        self._site_arrays = site_arrays
2366
+        lengths = [0] + [len(arr.tags) for arr in site_arrays]
2367
+        self._start_idxs = np.cumsum(lengths)[:-1]
2368
+
2369
+    def __len__(self):
2370
+        return sum(len(arr.tags) for arr in self._site_arrays)
2371
+
2372
+    def __getitem__(self, site):
2373
+        # treat tags as 1D sequence by defining custom dtype
2374
+        tag_dtype = np.asarray(site.tag).dtype
2375
+        dtype = np.dtype([('f{}'.format(i), tag_dtype)
2376
+                          for i in range(len(site.tag))])
2377
+        site_tag = np.asarray(site.tag).view(dtype)[0]
2378
+        # This slightly convoluted logic is necessary because there
2379
+        # may be >1 site_array with the same family for InfiniteSystems.
2380
+        for start, arr in zip(self._start_idxs, self._site_arrays):
2381
+            if arr.family != site.family:
2382
+                continue
2383
+            tags = arr.tags.view(dtype).reshape(-1)
2384
+            idx_in_fam = np.recarray.searchsorted(tags, site_tag)
2385
+            if idx_in_fam >= len(tags) or tags[idx_in_fam] != site_tag:
2386
+                continue
2387
+            return start + idx_in_fam
2388
+
2389
+        raise KeyError(site)
2390
+
2391
+
2392
+def _normalize_term_value(value, expected_n_values):
2393
+    try:
2394
+        value = np.asarray(value, dtype=complex)
2395
+    except TypeError:
2396
+        raise ValueError(
2397
+            "Matrix elements declared with incompatible shapes."
2398
+        ) from None
2399
+    # Upgrade to vector of matrices
2400
+    if len(value.shape) == 1:
2401
+        value = value[:, np.newaxis, np.newaxis]
2402
+    if len(value.shape) != 3:
2403
+        msg = (
2404
+            "Vectorized value functions must return an array of"
2405
+            "scalars or an array of matrices."
2406
+        )
2407
+        raise ValueError(msg)
2408
+    if value.shape[0] != expected_n_values:
2409
+        raise ValueError("Value functions must return a single value per "
2410
+                         "onsite/hopping.")
2411
+    return value
2412
+
2413
+
2414
+def _sort_term(term, value):
2415
+    term = np.asarray(term)
2416
+
2417
+    if not callable(value):
2418
+        value = _normalize_term_value(value, len(term))
2419
+        # Ensure that values still correspond to the correct
2420
+        # sites in 'term' once the latter has been sorted.
2421
+        value = value[term.argsort()]
2422
+
2423
+    term.sort()
2424
+
2425
+    return term, value
2426
+
2427
+
2428
+def _sort_hopping_term(term, value):
2429
+    term = np.asarray(term)
2430
+    # We want to sort the hopping terms in the same way
2431
+    # as tuples (i, j), so we need to use a record array.
2432
+    dtype = np.dtype([('f0', term.dtype), ('f1', term.dtype)])
2433
+    term_prime = term.view(dtype=dtype).reshape(-1)
2434
+    # _normalize_term will sort 'term' in-place via 'term_prime'
2435
+    _, value = _sort_term(term_prime, value)
2436
+
2437
+    return term, value
2438
+
2439
+
2440
+def _make_onsite_terms(builder, sites, site_offsets, term_offset):
2441
+    # Construct onsite terms.
2442
+    #
2443
+    # onsite_subgraphs
2444
+    #   Will contain tuples of the form described in
2445
+    #   kwant.system.VectorizedSystem.subgraphs, but during construction
2446
+    #   contains lists of the sites associated with each onsite term.
2447
+    #
2448
+    # onsite_term_values
2449
+    #   Contains a callable or array of constant values for each term.
2450
+    #
2451
+    # onsite_terms
2452
+    #   Constructed after the main loop. Contains a kwant.system.Term
2453
+    #   tuple for each onsite term.
2454
+    #
2455
+    # _onsite_term_by_site_id
2456
+    #   Maps the site ID to the number of the term that the site is
2457
+    #   a part of.
2458
+    #   lists the number of the
2459
+    #   Hamiltonian term associated with each site/hopping. For
2460
+    #   Hermitian conjugate hoppings "-term - 1" is stored instead.
2461
+
2462
+    onsite_subgraphs = []
2463
+    onsite_term_values = []
2464
+    onsite_term_parameters = []
2465
+    # We just use the cache to handle non/callables and exceptions
2466
+    cache = _value_params_pair_cache(1)
2467
+    # maps onsite key -> onsite term number
2468
+    onsite_to_term_nr = collections.OrderedDict()
2469
+    _onsite_term_by_site_id = []
2470
+    for site_id, site in enumerate(sites):
2471
+        val = builder.H[builder.symmetry.to_fd(site)][1]
2472
+        const_val = not callable(val)
2473
+        which_site_array = bisect.bisect(site_offsets, site_id) - 1
2474
+        # "key" uniquely identifies an onsite term.
2475
+        if const_val:
2476
+            key = (None, which_site_array)
2477
+        else:
2478
+            key = (id(val), which_site_array)
2479
+        if key not in onsite_to_term_nr:
2480
+            # Start a new term
2481
+            onsite_to_term_nr[key] = len(onsite_subgraphs)
2482
+            onsite_subgraphs.append([])
2483
+            if const_val:
2484
+                onsite_term_values.append([])
2485
+                onsite_term_parameters.append(None)
2486
+            else:
2487
+                val, params = cache(val)
2488
+                onsite_term_values.append(val)
2489
+                onsite_term_parameters.append(params)
2490
+        onsite_subgraphs[onsite_to_term_nr[key]].append(site_id)
2491
+        _onsite_term_by_site_id.append(onsite_to_term_nr[key])
2492
+        if const_val:
2493
+            vals = onsite_term_values[onsite_to_term_nr[key]]
2494
+            vals.append(val)
2495
+    # Sort the sites in each term, and normalize any constant
2496
+    # values to arrays of the correct dtype and shape.
2497
+    onsite_subgraphs, onsite_term_values = zip(*(
2498
+        _sort_term(term, value)
2499
+        for term, value in
2500
+        zip(onsite_subgraphs, onsite_term_values)))
2501
+    # Store site array index and site offsets rather than sites themselves
2502
+    tmp = []
2503
+    for (_, which), s in zip(onsite_to_term_nr, onsite_subgraphs):
2504
+        s = s - site_offsets[which]
2505
+        tmp.append(((which, which), (s, s)))
2506
+    onsite_subgraphs = tmp
2507
+    # onsite_term_errors[i] contains an exception if the corresponding
2508
+    # term has a value function with an illegal signature. We only raise
2509
+    # the exception if we actually try to evaluate the offending term
2510
+    # (to maintain backwards compatibility).
2511
+    onsite_term_errors = [
2512
+        err if isinstance(err, Exception) else None
2513
+        for err in onsite_term_parameters
2514
+    ]
2515
+    onsite_terms = [
2516
+        system.Term(
2517
+            subgraph=term_offset + i,
2518
+            hermitian=False,
2519
+            parameters=(
2520
+                params if not isinstance(params, Exception) else None
2521
+            ),
2522
+        )
2523
+        for i, params in enumerate(onsite_term_parameters)
2524
+    ]
2525
+    _onsite_term_by_site_id = np.array(_onsite_term_by_site_id)
2526
+
2527
+    return (onsite_subgraphs, onsite_terms, onsite_term_values,
2528
+            onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id)
2529
+
2530
+
2531
+def _make_hopping_terms(builder, graph, sites, site_offsets, cell_size, term_offset):
2532
+    # Construct the hopping terms
2533
+    #
2534
+    # The logic is the same as for the onsite terms, with the following
2535
+    # differences.
2536
+    #
2537
+    # _hopping_term_by_edge_id
2538
+    #   Maps hopping edge IDs to the number of the term that the hopping
2539
+    #   is a part of. For Hermitian conjugate hoppings "-term_number -1"
2540
+    #   is stored instead.
2541
+
2542
+    hopping_subgraphs = []
2543
+    hopping_term_values = []
2544
+    hopping_term_parameters = []
2545
+    # We just use the cache to handle non/callables and exceptions
2546
+    cache = _value_params_pair_cache(2)
2547
+    # maps hopping key -> hopping term number
2548
+    hopping_to_term_nr = collections.OrderedDict()
2549
+    _hopping_term_by_edge_id = []
2550
+    for tail, head in graph:
2551
+        tail_site, head_site = sites[tail], sites[head]
2552
+        if tail >= cell_size:
2553
+            # The tail belongs to the previous domain.  Find the
2554
+            # corresponding hopping with the tail in the fund. domain.
2555
+            tail_site, head_site = builder.symmetry.to_fd(tail_site, head_site)
2556
+        val = builder._get_edge(tail_site, head_site)
2557
+        herm_conj = val is Other
2558
+        if herm_conj:
2559
+            val = builder._get_edge(*builder.symmetry.to_fd(head_site, tail_site))
2560
+        const_val = not callable(val)
2561
+
2562
+        tail_site_array = bisect.bisect(site_offsets, tail) - 1
2563
+        head_site_array = bisect.bisect(site_offsets, head) - 1
2564
+        # "key" uniquely identifies a hopping term.
2565
+        if const_val:
2566
+            key = (None, tail_site_array, head_site_array)
2567
+        else:
2568
+            key = (id(val), tail_site_array, head_site_array)
2569
+        if herm_conj:
2570
+            key = (key[0], head_site_array, tail_site_array)
2571
+
2572
+        if key not in hopping_to_term_nr:
2573
+            # start a new term
2574
+            hopping_to_term_nr[key] = len(hopping_subgraphs)
2575
+            hopping_subgraphs.append([])
2576
+            if const_val:
2577
+                hopping_term_values.append([])
2578
+                hopping_term_parameters.append(None)
2579
+            else:
2580
+                val, params = cache(val)
2581
+                hopping_term_values.append(val)
2582
+                hopping_term_parameters.append(params)
2583
+
2584
+        if herm_conj:
2585
+            # Hopping terms come after all onsite terms, so we need
2586
+            # to offset the term ID by that many
2587
+            term_id = -term_offset - hopping_to_term_nr[key] - 1
2588
+            _hopping_term_by_edge_id.append(term_id)
2589
+        else:
2590
+            # Hopping terms come after all onsite terms, so we need
2591
+            # to offset the term ID by that many
2592
+            term_id = term_offset + hopping_to_term_nr[key]
2593
+            _hopping_term_by_edge_id.append(term_id)
2594
+            hopping_subgraphs[hopping_to_term_nr[key]].append((tail, head))
2595
+            if const_val:
2596
+                vals = hopping_term_values[hopping_to_term_nr[key]]
2597
+                vals.append(val)
2598
+    # Sort the hoppings in each term, and normalize any constant
2599
+    # values to arrays of the correct dtype and shape.
2600
+    if hopping_subgraphs:
2601
+        hopping_subgraphs, hopping_term_values = zip(*(
2602
+            _sort_hopping_term(term, value)
2603
+            for term, value in
2604
+            zip(hopping_subgraphs, hopping_term_values)))
2605
+    # Store site array index and site offsets rather than sites themselves
2606
+    tmp = []
2607
+    for (_, tail_which, head_which), h in zip(hopping_to_term_nr,
2608
+                                              hopping_subgraphs):
2609
+        start = (site_offsets[tail_which], site_offsets[head_which])
2610
+        # Transpose to get a pair of arrays rather than array of pairs
2611
+        # We use the fact that the underlying array is stored in
2612
+        # array-of-pairs order to search through it in 'hamiltonian'.
2613
+        tmp.append(((tail_which, head_which), (h - start).transpose()))
2614
+    hopping_subgraphs = tmp
2615
+    # hopping_term_errors[i] contains an exception if the corresponding
2616
+    # term has a value function with an illegal signature. We only raise
2617
+    # the exception if this becomes a problem (to maintain backwards
2618
+    # compatibility)
2619
+    hopping_term_errors = [
2620
+        err if isinstance(err, Exception) else None
2621
+        for err in hopping_term_parameters
2622
+    ]
2623
+    hopping_terms = [
2624
+        system.Term(
2625
+            subgraph=term_offset + i,
2626
+            hermitian=True,  # Builders are always Hermitian
2627
+            parameters=(
2628
+                params if not isinstance(params, Exception) else None
2629
+            ),
2630
+        )
2631
+        for i, params in enumerate(hopping_term_parameters)
2632
+    ]
2633
+    _hopping_term_by_edge_id = np.array(_hopping_term_by_edge_id)
2634
+
2635
+    return (hopping_subgraphs, hopping_terms, hopping_term_values,
2636
+            hopping_term_parameters, hopping_term_errors,
2637
+            _hopping_term_by_edge_id)
2638
+
2639
+
2640
+class FiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.FiniteVectorizedSystem):
2641
+    """Finalized `Builder` with leads.
2642
+
2643
+    Usable as input for the solvers in `kwant.solvers`.
2644
+
2645
+    Attributes
2646
+    ----------
2647
+    site_arrays : sequence of `~kwant.system.SiteArray`
2648
+        The sites of the system.
2649
+    sites : sequence
2650
+        ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2651
+        to the integer-labeled site ``i`` of the low-level system. The sites
2652
+        are ordered first by their family and then by their tag.
2653
+    id_by_site : mapping
2654
+        The inverse of ``sites``; maps high-level `~kwant.system.Site`
2655
+        instances to their integer label.
2656
+        Satisfies ``id_by_site[sites[i]] == i``.
2657
+    """
2658
+
2659
+    def __init__(self, builder):
2660
+        assert builder.symmetry.num_directions == 0
2661
+        assert builder.vectorize
2662
+
2663
+        sites = sorted(builder.H)
2664
+        id_by_site = {site: site_id for site_id, site in enumerate(sites)}
2665
+
2666
+        if not all(s.family.norbs for s in sites):
2667
+            raise ValueError("All sites must belong to families with norbs "
2668
+                             "specified for vectorized Builders. Specify "
2669
+                             "norbs when creating site families.")
2670
+
2671
+        graph = _make_graph(builder.H, id_by_site)
2672
+
2673
+        finalized_leads, lead_interfaces, lead_paddings =\
2674
+            _finalize_leads(builder.leads, id_by_site)
2675
+
2676
+        del id_by_site  # cleanup due to large size
2677
+
2678
+        site_arrays = _make_site_arrays(builder.H)
2679
+        # We need this to efficiently find which array a given
2680
+        # site belongs to
2681
+        site_offsets = np.cumsum([0] + [len(arr) for arr in site_arrays])
2682
+
2683
+        (onsite_subgraphs, onsite_terms, onsite_term_values,
2684
+         onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id) =\
2685
+            _make_onsite_terms(builder, sites, site_offsets, term_offset=0)
2686
+
2687
+        (hopping_subgraphs, hopping_terms, hopping_term_values,
2688
+         hopping_term_parameters, hopping_term_errors,
2689
+         _hopping_term_by_edge_id) =\
2690
+            _make_hopping_terms(builder, graph, sites, site_offsets,
2691
+                                len(sites), term_offset=len(onsite_terms))
2692
+
2693
+        # Construct the combined onsite/hopping term datastructures
2694
+        subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
2695
+        terms = tuple(onsite_terms) + tuple(hopping_terms)
2696
+        term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
2697
+        term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
2698
+
2699
+        # Construct system parameters
2700
+        parameters = set()
2701
+        for params in chain(onsite_term_parameters, hopping_term_parameters):
2702
+            if params is not None:
2703
+                parameters.update(params)
2704
+        parameters = frozenset(parameters)
2705
+
2706
+        self.site_arrays = site_arrays
2707
+        self.sites = _Sites(self.site_arrays)
2708
+        self.id_by_site = _IdBySite(self.site_arrays)
2709
+        self.graph = graph
2710
+        self.subgraphs = subgraphs
2711
+        self.terms = terms
2712
+        self._term_values = term_values
2713
+        self._term_errors = term_errors
2714
+        self._hopping_term_by_edge_id = _hopping_term_by_edge_id
2715
+        self._onsite_term_by_site_id = _onsite_term_by_site_id
2716
+        self.parameters = parameters
2717
+        self.symmetry = builder.symmetry
2718
+        self.leads = finalized_leads
2719
+        self.lead_interfaces = lead_interfaces
2720
+        self.lead_paddings = lead_paddings
2721
+        self._init_discrete_symmetries(builder)
2722
+
2723
+
2724
+def _make_lead_sites(builder, interface_order):
2725
+    #### For each site of the fundamental domain, determine whether it has
2726
+    #### neighbors in the previous domain or not.
2727
+    sym = builder.symmetry
2728
+    lsites_with = []       # Fund. domain sites with neighbors in prev. dom
2729
+    lsites_without = []    # Remaining sites of the fundamental domain
2730
+    for tail in builder.H: # Loop over all sites of the fund. domain.
2731
+        for head in builder._out_neighbors(tail):
2732
+            fd = sym.which(head)[0]
2733
+            if fd == 1:
2734
+                # Tail belongs to fund. domain, head to the next domain.
2735
+                lsites_with.append(tail)
2736
+                break
2737
+        else:
2738
+            # Tail is a fund. domain site not connected to prev. domain.
2739
+            lsites_without.append(tail)
2740
+
2741
+    if not lsites_with:
2742
+        warnings.warn('Infinite system with disconnected cells.',
2743
+                      RuntimeWarning, stacklevel=3)
2744
+
2745
+    ### Create list of sites and a lookup table
2746
+    minus_one = ta.array((-1,))
2747
+    plus_one = ta.array((1,))
2748
+    if interface_order is None:
2749
+        # interface must be sorted
2750
+        interface = [sym.act(minus_one, s) for s in lsites_with]
2751
+        interface.sort()
2752
+    else:
2753
+        lsites_with_set = set(lsites_with)
2754
+        lsites_with = []
2755
+        interface = []
2756
+        if interface_order:
2757
+            shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
2758
+        for shifted_iface_site in interface_order:
2759
+            # Shift the interface domain before the fundamental domain.
2760
+            # That's the right place for the interface of a lead to be, but
2761
+            # the sites of interface_order might live in a different
2762
+            # domain.
2763
+            iface_site = sym.act(shift, shifted_iface_site)
2764
+            lsite = sym.act(plus_one, iface_site)
2765
+
2766
+            try:
2767
+                lsites_with_set.remove(lsite)
2768
+            except KeyError:
2769
+                if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
2770
+                    raise ValueError(
2771
+                        'The sites in interface_order do not all '
2772
+                        'belong to the same lead cell.')
2773
+                else:
2774
+                    raise ValueError('A site in interface_order is not an '
2775
+                                     'interface site:\n' + str(iface_site))
2776
+            interface.append(iface_site)
2777
+            lsites_with.append(lsite)
2778
+        if lsites_with_set:
2779
+            raise ValueError(
2780
+                'interface_order did not contain all interface sites.')
2781
+        # `interface_order` *must* be sorted, hence `interface` should also
2782
+        if interface != sorted(interface):
2783
+            raise ValueError('Interface sites must be sorted.')
2784
+        del lsites_with_set
2785
+
2786
+    return sorted(lsites_with), sorted(lsites_without), interface
2787
+
2788
+
2789
+def _make_lead_graph(builder, sites, id_by_site, cell_size):
2790
+    sym = builder.symmetry
2791
+    g = graph.Graph()
2792
+    g.num_nodes = len(sites)  # Some sites could not appear in any edge.
2793
+    for tail_id, tail in enumerate(sites[:cell_size]):
2794
+        for head in builder._out_neighbors(tail):
2795
+            head_id = id_by_site.get(head)
2796
+            if head_id is None:
2797
+                # Head belongs neither to the fundamental domain nor to the
2798
+                # previous domain.  Check that it belongs to the next
2799
+                # domain and ignore it otherwise as an edge corresponding
2800
+                # to this one has been added already or will be added.
2801
+                fd = sym.which(head)[0]
2802
+                if fd != 1:
2803
+                    msg = ('Further-than-nearest-neighbor cells '
2804
+                           'are connected by hopping\n{0}.')
2805
+                    raise ValueError(msg.format((tail, head)))
2806
+                continue
2807
+            if head_id >= cell_size:
2808
+                # Head belongs to previous domain.  The edge added here
2809
+                # correspond to one left out just above.
2810
+                g.add_edge(head_id, tail_id)
2811
+            g.add_edge(tail_id, head_id)
2812
+
2813
+    return g.compressed()
2270 2814
 
2271 2815
 
2272 2816
 class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
... ...
@@ -2275,10 +2819,10 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2275 2819
     Attributes
2276 2820
     ----------
2277 2821
     sites : sequence
2278
-        ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
2822
+        ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2279 2823
         to the integer-labeled site ``i`` of the low-level system.
2280
-    id_by_site : dict
2281
-        The inverse of ``sites``; maps high-level `~kwant.builder.Site`
2824
+    id_by_site : mapping
2825
+        The inverse of ``sites``; maps high-level `~kwant.system.Site`
2282 2826
         instances to their integer label.
2283 2827
         Satisfies ``id_by_site[sites[i]] == i``.
2284 2828
 
... ...
@@ -2302,69 +2846,12 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2302 2846
         sym = builder.symmetry
2303 2847
         assert sym.num_directions == 1
2304 2848
 
2305
-        #### For each site of the fundamental domain, determine whether it has
2306
-        #### neighbors in the previous domain or not.
2307
-        lsites_with = []       # Fund. domain sites with neighbors in prev. dom
2308
-        lsites_without = []    # Remaining sites of the fundamental domain
2309
-        for tail in builder.H: # Loop over all sites of the fund. domain.
2310
-            for head in builder._out_neighbors(tail):
2311
-                fd = sym.which(head)[0]
2312
-                if fd == 1:
2313
-                    # Tail belongs to fund. domain, head to the next domain.
2314
-                    lsites_with.append(tail)
2315
-                    break
2316
-            else:
2317
-                # Tail is a fund. domain site not connected to prev. domain.
2318
-                lsites_without.append(tail)
2849
+        lsites_with, lsites_without, interface =\
2850
+            _make_lead_sites(builder, interface_order)
2319 2851
         cell_size = len(lsites_with) + len(lsites_without)
2320 2852
 
2321
-        if not lsites_with:
2322
-            warnings.warn('Infinite system with disconnected cells.',
2323
-                          RuntimeWarning, stacklevel=3)
2324
-
2325
-        ### Create list of sites and a lookup table
2326
-        minus_one = ta.array((-1,))
2327
-        plus_one = ta.array((1,))
2328
-        if interface_order is None:
2329
-            # interface must be sorted
2330
-            interface = [sym.act(minus_one, s) for s in lsites_with]
2331
-            interface.sort()
2332
-        else:
2333
-            lsites_with_set = set(lsites_with)
2334
-            lsites_with = []
2335
-            interface = []
2336
-            if interface_order:
2337
-                shift = ta.array((-sym.which(interface_order[0])[0] - 1,))
2338
-            for shifted_iface_site in interface_order:
2339
-                # Shift the interface domain before the fundamental domain.
2340
-                # That's the right place for the interface of a lead to be, but
2341
-                # the sites of interface_order might live in a different
2342
-                # domain.
2343
-                iface_site = sym.act(shift, shifted_iface_site)
2344
-                lsite = sym.act(plus_one, iface_site)
2345
-
2346
-                try:
2347
-                    lsites_with_set.remove(lsite)
2348
-                except KeyError:
2349
-                    if (-sym.which(shifted_iface_site)[0] - 1,) != shift:
2350
-                        raise ValueError(
2351
-                            'The sites in interface_order do not all '
2352
-                            'belong to the same lead cell.')
2353
-                    else:
2354
-                        raise ValueError('A site in interface_order is not an '
2355
-                                         'interface site:\n' + str(iface_site))
2356
-                interface.append(iface_site)
2357
-                lsites_with.append(lsite)
2358
-            if lsites_with_set:
2359
-                raise ValueError(
2360
-                    'interface_order did not contain all interface sites.')
2361
-            # `interface_order` *must* be sorted, hence `interface` should also
2362
-            if interface != sorted(interface):
2363
-                raise ValueError('Interface sites must be sorted.')
2364
-            del lsites_with_set
2365
-
2366 2853
         # we previously sorted the interface, so don't sort it again
2367
-        sites = sorted(lsites_with) + sorted(lsites_without) + interface
2854
+        sites = lsites_with + lsites_without + interface
2368 2855
         del lsites_with
2369 2856
         del lsites_without
2370 2857
         del interface
... ...
@@ -2372,41 +2859,20 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2372 2859
         for site_id, site in enumerate(sites):
2373 2860
             id_by_site[site] = site_id
2374 2861
 
2862
+        graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
2863
+
2375 2864
         # In the following, because many onsites/hoppings share the same
2376 2865
         # (value, parameter) pairs, we keep them in 'cache' so that we only
2377 2866
         # store a given pair in memory *once*. This is like interning strings.
2378 2867
 
2379
-        #### Make graph and extract onsite Hamiltonians.
2868
+        #### Extract onsites
2380 2869
         cache = _value_params_pair_cache(1)
2381
-        g = graph.Graph()
2382
-        g.num_nodes = len(sites)  # Some sites could not appear in any edge.
2383
-        onsites = []
2384
-        for tail_id, tail in enumerate(sites[:cell_size]):
2385
-            onsites.append(cache(builder.H[tail][1]))
2386
-            for head in builder._out_neighbors(tail):
2387
-                head_id = id_by_site.get(head)
2388
-                if head_id is None:
2389
-                    # Head belongs neither to the fundamental domain nor to the
2390
-                    # previous domain.  Check that it belongs to the next
2391
-                    # domain and ignore it otherwise as an edge corresponding
2392
-                    # to this one has been added already or will be added.
2393
-                    fd = sym.which(head)[0]
2394
-                    if fd != 1:
2395
-                        msg = ('Further-than-nearest-neighbor cells '
2396
-                               'are connected by hopping\n{0}.')
2397
-                        raise ValueError(msg.format((tail, head)))
2398
-                    continue
2399
-                if head_id >= cell_size:
2400
-                    # Head belongs to previous domain.  The edge added here
2401
-                    # correspond to one left out just above.
2402
-                    g.add_edge(head_id, tail_id)
2403
-                g.add_edge(tail_id, head_id)
2404
-        g = g.compressed()
2870
+        onsites = [cache(builder.H[tail][1]) for tail in sites[:cell_size]]
2405 2871
 
2406 2872
         #### Extract hoppings.
2407 2873
         cache = _value_params_pair_cache(2)
2408 2874
         hoppings = []
2409
-        for tail_id, head_id in g:
2875
+        for tail_id, head_id in graph:
2410 2876
             tail = sites[tail_id]
2411 2877
             head = sites[head_id]
2412 2878
             if tail_id >= cell_size:
... ...
@@ -2433,7 +2899,7 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2433 2899
         else:
2434 2900
             parameters = frozenset(parameters)
2435 2901
 
2436
-        self.graph = g
2902
+        self.graph = graph
2437 2903
         self.sites = sites
2438 2904
         self.site_ranges = _site_ranges(sites)
2439 2905
         self.id_by_site = id_by_site
... ...
@@ -2444,6 +2910,125 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2444 2910
         self.cell_size = cell_size
2445 2911
         self._init_discrete_symmetries(builder)
2446 2912
 
2913
+    def hamiltonian(self, i, j, *args, params=None):
2914
+        cs = self.cell_size
2915
+        if i == j >= cs:
2916
+            i -= cs
2917
+            j -= cs
2918
+        return super().hamiltonian(i, j, *args, params=params)
2919
+
2920
+
2921
+class InfiniteVectorizedSystem(_VectorizedFinalizedBuilderMixin, system.InfiniteVectorizedSystem):
2922
+    """Finalized infinite system, extracted from a `Builder`.
2923
+
2924
+    Attributes
2925
+    ----------
2926
+    site_arrays : sequence of `~kwant.system.SiteArray`
2927
+        The sites of the system.
2928
+    sites : sequence
2929
+        ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
2930
+        to the integer-labeled site ``i`` of the low-level system.
2931
+    id_by_site : mapping
2932
+        The inverse of ``sites``; maps high-level `~kwant.system.Site`
2933
+        instances to their integer label.
2934
+        Satisfies ``id_by_site[sites[i]] == i``.
2935
+
2936
+    Notes
2937
+    -----
2938
+    In infinite systems ``sites`` and ``site_arrays`` consists of 3 parts:
2939
+    sites in the fundamental domain (FD) with hoppings to neighboring cells,
2940
+    sites in the FD with no hoppings to neighboring cells, and sites in FD+1
2941
+    attached to the FD by hoppings. Each of these three subsequences is
2942
+    individually sorted.
2943
+    """
2944
+
2945
+    def __init__(self, builder, interface_order=None):
2946
+        """
2947
+        Finalize a builder instance which has to have exactly a single
2948
+        symmetry direction.
2949
+
2950
+        If interface_order is not set, the order of the interface sites in the
2951
+        finalized system will be arbitrary.  If interface_order is set to a
2952
+        sequence of interface sites, this order will be kept.
2953
+        """
2954
+        sym = builder.symmetry
2955
+        assert sym.num_directions == 1
2956
+        assert builder.vectorize
2957
+
2958
+        lsites_with, lsites_without, interface =\
2959
+            _make_lead_sites(builder, interface_order)
2960
+        cell_size = len(lsites_with) + len(lsites_without)
2961
+
2962
+
2963
+        sites = lsites_with + lsites_without + interface
2964
+        id_by_site = {}
2965
+        for site_id, site in enumerate(sites):
2966
+            id_by_site[site] = site_id
2967
+        # these are needed for constructing hoppings
2968
+        lsites_with = frozenset(lsites_with)
2969
+        lsites_without = frozenset(lsites_without)
2970
+        interface = frozenset(interface)
2971
+
2972
+        if not all(s.family.norbs for s in sites):
2973
+            raise ValueError("All sites must belong to families with norbs "
2974
+                             "specified for vectorized Builders. Specify "
2975
+                             "norbs when creating site families.")
2976
+
2977
+        graph = _make_lead_graph(builder, sites, id_by_site, cell_size)
2978
+
2979
+        del id_by_site  # cleanup due to large size
2980
+
2981
+        # In order to conform to the kwant.system.InfiniteVectorizedSystem
2982
+        # interface we need to put the sites that connect to the previous
2983
+        # cell *first*, then the sites that do not couple to the previous
2984
+        # cell, then the sites in the previous cell. Because sites in
2985
+        # a SiteArray are sorted by tag this means that the sites in these
2986
+        # 3 different sets need to be in different SiteArrays.
2987
+        site_arrays = (
2988
+            _make_site_arrays(lsites_with)
2989
+            + _make_site_arrays(lsites_without)
2990
+            + _make_site_arrays(interface)
2991
+        )
2992
+
2993
+        site_offsets = np.cumsum([0] + [len(arr) for arr in site_arrays])
2994
+
2995
+        (onsite_subgraphs, onsite_terms, onsite_term_values,
2996
+         onsite_term_parameters, onsite_term_errors, _onsite_term_by_site_id) =\
2997
+            _make_onsite_terms(builder, sites, site_offsets, term_offset=0)
2998
+
2999
+        (hopping_subgraphs, hopping_terms, hopping_term_values,
3000
+         hopping_term_parameters, hopping_term_errors,
3001
+         _hopping_term_by_edge_id) =\
3002
+            _make_hopping_terms(builder, graph, sites, site_offsets,
3003
+                                cell_size, term_offset=len(onsite_terms))
3004
+
3005
+        # Construct the combined onsite/hopping term datastructures
3006
+        subgraphs = tuple(onsite_subgraphs) + tuple(hopping_subgraphs)
3007
+        terms = tuple(onsite_terms) + tuple(hopping_terms)
3008
+        term_values = tuple(onsite_term_values) + tuple(hopping_term_values)
3009
+        term_errors = tuple(onsite_term_errors) + tuple(hopping_term_errors)
3010
+
3011
+        # Construct system parameters
3012
+        parameters = set()
3013
+        for params in chain(onsite_term_parameters, hopping_term_parameters):
3014
+            if params is not None:
3015
+                parameters.update(params)
3016
+        parameters = frozenset(parameters)
3017
+
3018
+        self.site_arrays = site_arrays
3019
+        self.sites = _Sites(self.site_arrays)
3020
+        self.id_by_site = _IdBySite(self.site_arrays)
3021
+        self.graph = graph
3022
+        self.subgraphs = subgraphs
3023
+        self.terms = terms
3024
+        self._term_values = term_values
3025
+        self._term_errors = term_errors
3026
+        self._hopping_term_by_edge_id = _hopping_term_by_edge_id
3027
+        self._onsite_term_by_site_id = _onsite_term_by_site_id
3028
+        self.parameters = parameters
3029
+        self.symmetry = builder.symmetry
3030
+        self.cell_size = cell_size
3031
+        self._init_discrete_symmetries(builder)
2447 3032
 
2448 3033
     def hamiltonian(self, i, j, *args, params=None):
2449 3034
         cs = self.cell_size
... ...
@@ -2452,5 +3037,15 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
2452 3037
             j -= cs
2453 3038
         return super().hamiltonian(i, j, *args, params=params)
2454 3039
 
2455
-    def pos(self, i):
2456
-        return self.sites[i].pos
3040
+
3041
+def is_finite_system(syst):
3042
+    return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
3043
+
3044
+
3045
+def is_infinite_system(syst):
3046
+    return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
3047
+
3048
+
3049
+def is_system(syst):
3050
+    return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem,
3051
+                             InfiniteSystem, InfiniteVectorizedSystem))
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2017 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -20,7 +20,7 @@ from sympy.printing.lambdarepr import LambdaPrinter
20 20
 from sympy.printing.precedence import precedence
21 21
 from sympy.core.function import AppliedUndef
22 22
 
23
-from .. import builder, lattice
23
+from .. import builder, lattice, system
24 24
 from .. import KwantDeprecationWarning
25 25
 from .._common import reraise_warnings
26 26
 from ._common import (sympify, gcd, position_operators, momentum_operators,
... ...
@@ -63,7 +63,7 @@ class _DiscretizedBuilder(builder.Builder):
63 63
 
64 64
         for key, val in itertools.chain(self.site_value_pairs(),
65 65
                                         self.hopping_value_pairs()):
66
-            if isinstance(key, builder.Site):
66
+            if isinstance(key, system.Site):
67 67
                 result.append("# Onsite element:\n")
68 68
             else:
69 69
                 a, b = key
... ...
@@ -176,7 +176,7 @@ class LandauLattice(kwant.lattice.Monatomic):
176 176
     """
177 177
     A `~kwant.lattice.Monatomic` lattice with a Landau level index per site.
178 178
 
179
-    Site tags (see `~kwant.builder.SiteFamily`) are pairs of integers, where
179
+    Site tags (see `~kwant.system.SiteFamily`) are pairs of integers, where
180 180
     the first integer describes the real space position and the second the
181 181
     Landau level index.
182 182
 
... ...
@@ -1,5 +1,5 @@
1 1
 # -*- coding: utf-8 -*-
2
-# Copyright 2011-2016 Kwant authors.
2
+# Copyright 2011-2019 Kwant authors.
3 3
 #
4 4
 # This file is part of Kwant.  It is subject to the license terms in the file
5 5
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -928,7 +928,7 @@ def RandomVectors(syst, where=None, rng=None):
928 928
     ----------
929 929
     syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
930 930
         If a system is passed, it should contain no leads.
931
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
931
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
932 932
         Spatial range of the vectors produced. If ``syst`` is a
933 933
         `~kwant.system.FiniteSystem`, where behaves as in
934 934
         `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
... ...
@@ -950,7 +950,7 @@ class LocalVectors:
950 950
     ----------
951 951
     syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
952 952
         If a system is passed, it should contain no leads.
953
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
953
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
954 954
         Spatial range of the vectors produced. If ``syst`` is a
955 955
         `~kwant.system.FiniteSystem`, where behaves as in
956 956
         `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2013 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -13,7 +13,7 @@ from math import sqrt
13 13
 from itertools import product
14 14
 import numpy as np
15 15
 import tinyarray as ta
16
-from . import builder
16
+from . import builder, system
17 17
 from .linalg import lll
18 18
 from ._common import ensure_isinstance
19 19
 
... ...
@@ -70,7 +70,7 @@ class Polyatomic:
70 70
     A Bravais lattice with an arbitrary number of sites in the basis.
71 71
 
72 72
     Contains `Monatomic` sublattices.  Note that an instance of ``Polyatomic`` is
73
-    not itself a `~kwant.builder.SiteFamily`, only its sublattices are.
73
+    not itself a `~kwant.system.SiteFamily`, only its sublattices are.
74 74
 
75 75
     Parameters
76 76
     ----------
... ...
@@ -171,7 +171,7 @@ class Polyatomic:
171 171
         >>> syst[lat.neighbors()] = 1
172 172
         """
173 173
         def shape_sites(symmetry=None):
174
-            Site = builder.Site
174
+            Site = system.Site
175 175
 
176 176
             if symmetry is None:
177 177
                 symmetry = builder.NoSymmetry()
... ...
@@ -306,7 +306,7 @@ class Polyatomic:
306 306
         """
307 307
         # This algorithm is not designed to be fast and can be improved,
308 308
         # however there is no real need.
309
-        Site = builder.Site
309
+        Site = system.Site
310 310
         sls = self.sublattices
311 311
         shortest_hopping = sls[0].n_closest(
312 312
             sls[0].pos(([0] * sls[0].lattice_dim)), 2)[-1]
... ...
@@ -396,12 +396,12 @@ def short_array_str(array):
396 396
     return full[1 : -1]
397 397
 
398 398
 
399
-class Monatomic(builder.SiteFamily, Polyatomic):
399
+class Monatomic(system.SiteFamily, Polyatomic):
400 400
     """
401 401
     A Bravais lattice with a single site in the basis.
402 402
 
403
-    Instances of this class provide the `~kwant.builder.SiteFamily` interface.
404
-    Site tags (see `~kwant.builder.SiteFamily`) are sequences of integers and
403
+    Instances of this class provide the `~kwant.system.SiteFamily` interface.
404
+    Site tags (see `~kwant.system.SiteFamily`) are sequences of integers and
405 405
     describe the lattice coordinates of a site.
406 406
 
407 407
     ``Monatomic`` instances are used as site families on their own or as
... ...
@@ -470,6 +470,13 @@ class Monatomic(builder.SiteFamily, Polyatomic):
470 470
     def __str__(self):
471 471
         return self.cached_str
472 472
 
473
+    def normalize_tags(self, tags):
474
+        tags = np.asarray(tags, int)
475
+        tags.flags.writeable = False
476
+        if tags.shape[1] != self.lattice_dim:
477
+            raise ValueError("Dimensionality mismatch.")
478
+        return tags
479
+
473 480
     def normalize_tag(self, tag):
474 481
         tag = ta.array(tag, int)
475 482
         if len(tag) != self.lattice_dim:
... ...
@@ -495,6 +502,10 @@ class Monatomic(builder.SiteFamily, Polyatomic):
495 502
         """
496 503
         return ta.array(self.n_closest(pos)[0])
497 504
 
505
+    def positions(self, tags):
506
+        """Return the real-space positions of the sites with the given tags."""
507
+        return tags @ self._prim_vecs + self.offset
508
+
498 509
     def pos(self, tag):
499 510
         """Return the real-space position of the site with a given tag."""
500 511
         return ta.dot(tag, self._prim_vecs) + self.offset
... ...
@@ -687,26 +698,48 @@ class TranslationalSymmetry(builder.Symmetry):
687 698
 
688 699
     def which(self, site):
689 700
         det_x_inv_m_part, det_m = self._get_site_family_data(site.family)[-2:]
690
-        result = ta.dot(det_x_inv_m_part, site.tag) // det_m
701
+        if isinstance(site, system.Site):
702
+            result = ta.dot(det_x_inv_m_part, site.tag) // det_m
703
+        elif isinstance(site, system.SiteArray):
704
+            result = np.dot(det_x_inv_m_part, site.tags.transpose()) // det_m
705
+        else:
706
+            raise TypeError("'site' must be a Site or a SiteArray")
707
+
691 708
         return -result if self.is_reversed else result
692 709
 
693 710
     def act(self, element, a, b=None):
694
-        element = ta.array(element)
695
-        if element.dtype is not int:
711
+        is_site = isinstance(a, system.Site)
712
+        # Tinyarray for small arrays (single site) else numpy
713
+        array_mod = ta if is_site else np
714
+        element = array_mod.array(element)
715
+        if not np.issubdtype(element.dtype, np.integer):
696 716
             raise ValueError("group element must be a tuple of integers")
717
+        if (len(element.shape) == 2 and is_site):
718
+            raise ValueError("must provide a single group element when "
719
+                             "acting on single sites.")
720
+        if (len(element.shape) == 1 and not is_site):
721
+            raise ValueError("must provide a sequence of group elements "
722
+                             "when acting on site arrays.")
697 723
         m_part = self._get_site_family_data(a.family)[0]
698 724
         try:
699
-            delta = ta.dot(m_part, element)
725
+            delta = array_mod.dot(m_part, element)
700 726
         except ValueError:
701 727
             msg = 'Expecting a {0}-tuple group element, but got `{1}` instead.'
702 728
             raise ValueError(msg.format(self.num_directions, element))
703 729
         if self.is_reversed:
704 730
             delta = -delta
705 731
         if b is None:
706
-            return builder.Site(a.family, a.tag + delta, True)
732
+            if is_site:
733
+                return system.Site(a.family, a.tag + delta, True)
734
+            else:
735
+                return system.SiteArray(a.family, a.tags + delta.transpose())
707 736
         elif b.family == a.family:
708
-            return (builder.Site(a.family, a.tag + delta, True),
709
-                    builder.Site(b.family, b.tag + delta, True))
737
+            if is_site:
738
+                return (system.Site(a.family, a.tag + delta, True),
739
+                        system.Site(b.family, b.tag + delta, True))
740
+            else:
741
+                return (system.SiteArray(a.family, a.tags + delta.transpose()),
742
+                        system.SiteArray(b.family, b.tags + delta.transpose()))
710 743
         else:
711 744
             m_part = self._get_site_family_data(b.family)[0]
712 745
             try:
... ...
@@ -717,8 +750,12 @@ class TranslationalSymmetry(builder.Symmetry):
717 750
                 raise ValueError(msg.format(self.num_directions, element))
718 751
             if self.is_reversed:
719 752
                 delta2 = -delta2
720
-            return (builder.Site(a.family, a.tag + delta, True),
721
-                    builder.Site(b.family, b.tag + delta2, True))
753
+            if is_site:
754
+                return (system.Site(a.family, a.tag + delta, True),
755
+                        system.Site(b.family, b.tag + delta2, True))
756
+            else:
757
+                return (system.SiteArray(a.family, a.tags + delta.transpose()),
758
+                        system.SiteArray(b.family, b.tags + delta2.transpose()))
722 759
 
723 760
     def reversed(self):
724 761
         """Return a reversed copy of the symmetry.
... ...
@@ -13,7 +13,7 @@ import cython
13 13
 from operator import itemgetter
14 14
 import functools as ft
15 15
 import collections
16
-import numbers
16
+import warnings
17 17
 
18 18
 import numpy as np
19 19
 import tinyarray as ta
... ...
@@ -25,8 +25,10 @@ from .graph.core cimport EdgeIterator
25 25
 from .graph.core import DisabledFeatureError, NodeDoesNotExistError
26 26
 from .graph.defs cimport gint
27 27
 from .graph.defs import gint_dtype
28
-from .system import InfiniteSystem
29
-from ._common import UserCodeError, get_parameters, deprecate_args
28
+from .system import is_infinite, Site
29
+from ._common import (
30
+    UserCodeError, KwantDeprecationWarning, get_parameters, deprecate_args
31
+)
30 32
 
31 33
 
32 34
 ################ Generic Utility functions
... ...
@@ -150,7 +152,7 @@ def _get_all_orbs(gint[:, :] where, gint[:, :] site_ranges):
150 152
 
151 153
 def _get_tot_norbs(syst):
152 154
     cdef gint _unused, tot_norbs
153
-    is_infinite_system = isinstance(syst, InfiniteSystem)
155
+    is_infinite_system = is_infinite(syst)
154 156
     n_sites = syst.cell_size if is_infinite_system else syst.graph.num_nodes
155 157
     _get_orbs(np.asarray(syst.site_ranges, dtype=gint_dtype),
156 158
               n_sites, &tot_norbs, &_unused)
... ...
@@ -166,7 +168,7 @@ def _normalize_site_where(syst, where):
166 168
     otherwise it should contain integers.
167 169
     """
168 170
     if where is None:
169
-        if isinstance(syst, InfiniteSystem):
171
+        if is_infinite(syst):
170 172
             where = list(range(syst.cell_size))
171 173
         else:
172 174
             where = list(range(syst.graph.num_nodes))
... ...
@@ -174,13 +176,12 @@ def _normalize_site_where(syst, where):
174 176
         try:
175 177
             where = [syst.id_by_site[s] for s in filter(where, syst.sites)]
176 178
         except AttributeError:
177
-            if isinstance(syst, InfiniteSystem):
179
+            if is_infinite(syst):
178 180
                 where = [s for s in range(syst.cell_size) if where(s)]
179 181
             else:
180 182
                 where = [s for s in range(syst.graph.num_nodes) if where(s)]
181 183
     else:
182
-        # Cannot check for builder.Site due to circular imports
183
-        if not isinstance(where[0], numbers.Integral):
184
+        if isinstance(where[0], Site):
184 185
             try:
185 186
                 where = [syst.id_by_site[s] for s in where]
186 187
             except AttributeError:
... ...
@@ -189,7 +190,7 @@ def _normalize_site_where(syst, where):
189 190
 
190 191
     where = np.asarray(where, dtype=gint_dtype).reshape(-1, 1)
191 192
 
192
-    if isinstance(syst, InfiniteSystem) and np.any(where >= syst.cell_size):
193
+    if is_infinite(syst) and np.any(where >= syst.cell_size):
193 194
         raise ValueError('Only sites in the fundamental domain may be '
194 195
                          'specified using `where`.')
195 196
     if np.any(np.logical_or(where < 0, where >= syst.graph.num_nodes)):
... ...
@@ -210,7 +211,7 @@ def _normalize_hopping_where(syst, where):
210 211
     if where is None:
211 212
         # we cannot extract the hoppings in the same order as they are in the
212 213
         # graph while simultaneously excluding all inter-cell hoppings
213
-        if isinstance(syst, InfiniteSystem):
214
+        if is_infinite(syst):
214 215
             raise ValueError('`where` must be provided when calculating '
215 216
                              'current in an InfiniteSystem.')
216 217
         where = list(syst.graph)
... ...
@@ -223,8 +224,7 @@ def _normalize_hopping_where(syst, where):
223 224
         else:
224 225
             where = list(filter(lambda h: where(*h), syst.graph))
225 226
     else:
226
-        # Cannot check for builder.Site due to circular imports
227
-        if not isinstance(where[0][0], numbers.Integral):
227
+        if isinstance(where[0][0], Site):
228 228
             try:
229 229
                 where = list((syst.id_by_site[a], syst.id_by_site[b])
230 230
                                for a, b in where)
... ...
@@ -244,7 +244,7 @@ def _normalize_hopping_where(syst, where):
244 244
 
245 245
     where = np.asarray(where, dtype=gint_dtype)
246 246
 
247
-    if isinstance(syst, InfiniteSystem) and np.any(where > syst.cell_size):
247
+    if is_infinite(syst) and np.any(where > syst.cell_size):
248 248
         raise ValueError('Only intra-cell hoppings may be specified '
249 249
                          'using `where`.')
250 250
 
... ...
@@ -702,7 +702,11 @@ cdef class _LocalOperator:
702 702
             return mat
703 703
 
704 704
         offsets, norbs = _get_all_orbs(self.where, self._site_ranges)
705
-        return  BlockSparseMatrix(self.where, offsets, norbs, get_ham)
705
+        # TODO: update operators to use 'hamiltonian_term' rather than
706
+        #       'hamiltonian'.
707
+        with warnings.catch_warnings():
708
+            warnings.simplefilter("ignore", category=KwantDeprecationWarning)
709
+            return  BlockSparseMatrix(self.where, offsets, norbs, get_ham)
706 710
 
707 711
     def __getstate__(self):
708 712
         return (
... ...
@@ -735,10 +739,10 @@ cdef class Density(_LocalOperator):
735 739
         maps from site families to square matrices. If a function is given it
736 740
         must take the same arguments as the onsite Hamiltonian functions of the
737 741
         system.
738
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
742
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
739 743
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
740 744
         then this should be a sequence of integers. If a function is provided,
741
-        it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
745
+        it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
742 746
         a finalized builder) and return True or False.  If not provided, the
743 747
         operator will be calculated over all sites in the system.
744 748
     check_hermiticity: bool
... ...
@@ -884,11 +888,11 @@ cdef class Current(_LocalOperator):
884 888
         matrices (scalars are allowed if the site family has 1 orbital per
885 889
         site). If a function is given it must take the same arguments as the
886 890
         onsite Hamiltonian functions of the system.
887
-    where : sequence of pairs of `int` or `~kwant.builder.Site`, or callable, optional
891
+    where : sequence of pairs of `int` or `~kwant.system.Site`, or callable, optional
888 892
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
889 893
         then this should be a sequence of pairs of integers. If a function is
890 894
         provided, it should take a pair of integers or a pair of
891
-        `~kwant.builder.Site` (if ``syst`` is a finalized builder) and return
895
+        `~kwant.system.Site` (if ``syst`` is a finalized builder) and return
892 896
         True or False.  If not provided, the operator will be calculated over
893 897
         all hoppings in the system.
894 898
     check_hermiticity : bool
... ...
@@ -1016,10 +1020,10 @@ cdef class Source(_LocalOperator):
1016 1020
         matrices (scalars are allowed if the site family has 1 orbital per
1017 1021
         site). If a function is given it must take the same arguments as the
1018 1022
         onsite Hamiltonian functions of the system.
1019
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
1023
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
1020 1024
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
1021 1025
         then this should be a sequence of integers. If a function is provided,
1022
-        it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
1026
+        it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
1023 1027
         a finalized builder) and return True or False.  If not provided, the
1024 1028
         operator will be calculated over all sites in the system.
1025 1029
     check_hermiticity : bool
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2018 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -303,7 +303,7 @@ def _loops_in_infinite(syst):
303 303
     extended_sites : callable : int -> Site
304 304
         Given a site index in the extended system consisting of
305 305
         two unit cells, returns the associated high-level
306
-        `kwant.builder.Site`.
306
+        `kwant.system.Site`.
307 307
     """
308 308
     assert isinstance(syst, system.InfiniteSystem)
309 309
     _check_infinite_syst(syst)
... ...
@@ -375,7 +375,7 @@ def _loops_in_composite(syst):
375 375
         -1 if the site is part of the reduced scattering region (see notes).
376 376
     extended_sites : callable : int -> Site
377 377
         Given a site index in the extended scattering region (see notes),
378
-        returns the associated high-level `kwant.builder.Site`.
378
+        returns the associated high-level `kwant.system.Site`.
379 379
 
380 380
     Notes
381 381
     -----
... ...
@@ -401,7 +401,7 @@ def _loops_in_composite(syst):
401 401
     # Get distance matrix for the extended scattering region,
402 402
     # a function that maps sites to their lead patches (-1 for sites
403 403
     # in the reduced scattering region), and a function that maps sites
404
-    # to high-level 'kwant.builder.Site' objects.
404
+    # to high-level 'kwant.system.Site' objects.
405 405
     distance_matrix, which_patch, extended_sites =\
406 406
         _extended_scattering_region(syst)
407 407
 
... ...
@@ -439,7 +439,7 @@ def _extended_scattering_region(syst):
439 439
         -1 if the site is part of the reduced scattering region.
440 440
     extended_sites : callable : int -> Site
441 441
         Given a site index in the extended scattering region, returns
442
-        the associated high-level `kwant.builder.Site`.
442
+        the associated high-level `kwant.system.Site`.
443 443
 
444 444
     Notes
445 445
     -----
... ...
@@ -1027,7 +1027,7 @@ class magnetic_gauge:
1027 1027
             The first callable computes the Peierls phase in the scattering
1028 1028
             region and the remaining callables compute the Peierls phases
1029 1029
             in each of the leads. Each callable takes a pair of
1030
-            `~kwant.builder.Site` (a hopping) and returns a unit complex
1030
+            `~kwant.system.Site` (a hopping) and returns a unit complex
1031 1031
             number (Peierls phase) that multiplies that hopping.
1032 1032
         """
1033 1033
         return self._peierls(syst_field, *lead_fields, tol=tol, average=False)
... ...
@@ -1,5 +1,5 @@
1 1
 # -*- coding: utf-8 -*-
2
-# Copyright 2011-2018 Kwant authors.
2
+# Copyright 2011-2019 Kwant authors.
3 3
 #
4 4
 # This file is part of Kwant.  It is subject to the license terms in the file
5 5
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -402,7 +402,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
402 402
     Returns
403 403
     -------
404 404
     sites : list of (site, lead_number, copy_number) tuples
405
-        A site is a `~kwant.builder.Site` instance if the system is not finalized,
405
+        A site is a `~kwant.system.Site` instance if the system is not finalized,
406 406
         and an integer otherwise.  For system sites `lead_number` is `None` and
407 407
         `copy_number` is `0`, for leads both are integers.
408 408
     lead_cells : list of slices
... ...
@@ -427,7 +427,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
427 427
                               lead.builder.sites() for i in
428 428
                               range(num_lead_cells)))
429 429
             lead_cells.append(slice(start, len(sites)))
430
-    elif isinstance(syst, system.FiniteSystem):
430
+    elif system.is_finite(syst):
431 431
         sites = [(i, None, 0) for i in range(syst.graph.num_nodes)]
432 432
         for leadnr, lead in enumerate(syst.leads):
433 433
             start = len(sites)
... ...
@@ -528,7 +528,7 @@ def sys_leads_hoppings(sys, num_lead_cells=2):
528 528
     Returns
529 529
     -------
530 530
     hoppings : list of (hopping, lead_number, copy_number) tuples
531
-        A site is a `~kwant.builder.Site` instance if the system is not finalized,
531
+        A site is a `~kwant.system.Site` instance if the system is not finalized,
532 532
         and an integer otherwise.  For system sites `lead_number` is `None` and
533 533
         `copy_number` is `0`, for leads both are integers.
534 534
     lead_cells : list of slices
... ...
@@ -812,7 +812,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
812 812
 
813 813
     - The meaning of "site" depends on whether the system to be plotted is a
814 814
       builder or a low level system.  For builders, a site is a
815
-      kwant.builder.Site object.  For low level systems, a site is an integer
815
+      kwant.system.Site object.  For low level systems, a site is an integer
816 816
       -- the site number.
817 817
 
818 818
     - color and symbol definitions may be tuples, but not lists or arrays.
... ...
@@ -961,7 +961,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
961 961
 
962 962
     if site_color is None:
963 963
         cycle = _color_cycle()
964
-        if isinstance(syst, (builder.FiniteSystem, builder.InfiniteSystem)):
964
+        if builder.is_system(syst):
965 965
             # Skipping the leads for brevity.
966 966
             families = sorted({site.family for site in syst.sites})
967 967
             color_mapping = dict(zip(families, cycle))
... ...
@@ -1291,7 +1291,7 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
1291 1291
     if callable(value):
1292 1292
         value = [value(site[0]) for site in sites]
1293 1293
     else:
1294
-        if not isinstance(syst, system.FiniteSystem):
1294
+        if not system.is_finite(syst):
1295 1295
             raise ValueError('List of values is only allowed as input '
1296 1296
                              'for finalized systems.')
1297 1297
     value = np.array(value)
... ...
@@ -1407,7 +1407,7 @@ def bands(sys, args=(), momenta=65, file=None, show=True, dpi=None,
1407 1407
                            "for bands()")
1408 1408
 
1409 1409
     syst = sys  # for naming consistency inside function bodies
1410
-    _common.ensure_isinstance(syst, system.InfiniteSystem)
1410
+    _common.ensure_isinstance(syst, (system.InfiniteSystem, system.InfiniteVectorizedSystem))
1411 1411
 
1412 1412
     momenta = np.array(momenta)
1413 1413
     if momenta.ndim != 1:
... ...
@@ -1483,7 +1483,7 @@ def spectrum(syst, x, y=None, params=None, mask=None, file=None,
1483 1483
     if y is not None and not _p.has3d:
1484 1484
         raise RuntimeError("Installed matplotlib does not support 3d plotting")
1485 1485
 
1486
-    if isinstance(syst, system.FiniteSystem):
1486
+    if system.is_finite(syst):
1487 1487
         def ham(**kwargs):
1488 1488
             return syst.hamiltonian_submatrix(params=kwargs, sparse=False)
1489 1489
     elif callable(syst):
... ...
@@ -1751,7 +1751,7 @@ def interpolate_current(syst, current, relwidth=None, abswidth=None, n=9):
1751 1751
         the extents of `field`: ((x0, x1), (y0, y1), ...)
1752 1752
 
1753 1753
     """
1754
-    if not isinstance(syst, builder.FiniteSystem):
1754
+    if not builder.is_finite_system(syst):
1755 1755
         raise TypeError("The system needs to be finalized.")
1756 1756
 
1757 1757
     if len(current) != syst.graph.num_edges:
... ...
@@ -1844,7 +1844,7 @@ def interpolate_density(syst, density, relwidth=None, abswidth=None, n=9,
1844 1844
         the extents of ``field``: ((x0, x1), (y0, y1), ...)
1845 1845
 
1846 1846
     """
1847
-    if not isinstance(syst, builder.FiniteSystem):
1847
+    if not builder.is_finite_system(syst):
1848 1848
         raise TypeError("The system needs to be finalized.")
1849 1849
 
1850 1850
     if len(density) != len(syst.sites):
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2013 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -8,13 +8,269 @@
8 8
 
9 9
 """Low-level interface of systems"""
10 10
 
11
-__all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
11
+__all__ = [
12
+    'Site', 'SiteArray', 'SiteFamily',
13
+    'System', 'VectorizedSystem', 'FiniteSystem', 'FiniteVectorizedSystem',
14
+    'InfiniteSystem', 'InfiniteVectorizedSystem',
15
+]
12 16
 
13 17
 import abc
14 18
 import warnings
19
+import operator
15 20
 from copy import copy
21
+from collections import namedtuple
22
+from functools import total_ordering, lru_cache
23
+import numpy as np
16 24
 from . import _system
17
-from ._common  import deprecate_args
25
+from ._common  import deprecate_args, KwantDeprecationWarning
26
+
27
+
28
+
29
+################ Sites and Site families
30
+
31
+class Site(tuple):
32
+    """A site, member of a `SiteFamily`.
33
+
34
+    Sites are the vertices of the graph which describes the tight binding
35
+    system in a `Builder`.
36
+
37
+    A site is uniquely identified by its family and its tag.
38
+
39
+    Parameters
40
+    ----------
41
+    family : an instance of `SiteFamily`
42
+        The 'type' of the site.
43
+    tag : a hashable python object
44
+        The unique identifier of the site within the site family, typically a
45
+        vector of integers.
46
+
47
+    Raises
48
+    ------
49
+    ValueError
50
+        If `tag` is not a proper tag for `family`.
51
+
52
+    Notes
53
+    -----
54
+    For convenience, ``family(*tag)`` can be used instead of ``Site(family,
55
+    tag)`` to create a site.
56
+
57
+    The parameters of the constructor (see above) are stored as instance
58
+    variables under the same names.  Given a site ``site``, common things to
59
+    query are thus ``site.family``, ``site.tag``, and ``site.pos``.
60
+    """
61
+    __slots__ = ()
62
+
63
+    family = property(operator.itemgetter(0),
64
+                      doc="The site family to which the site belongs.")
65
+    tag = property(operator.itemgetter(1), doc="The tag of the site.")
66
+
67
+
68
+    def __new__(cls, family, tag, _i_know_what_i_do=False):
69
+        if _i_know_what_i_do:
70
+            return tuple.__new__(cls, (family, tag))
71
+        try:
72
+            tag = family.normalize_tag(tag)
73
+        except (TypeError, ValueError) as e:
74
+            msg = 'Tag {0} is not allowed for site family {1}: {2}'
75
+            raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
76
+        return tuple.__new__(cls, (family, tag))
77
+
78
+    def __repr__(self):
79
+        return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
80
+
81
+    def __str__(self):
82
+        sf = self.family
83
+        return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
84
+
85
+    def __getnewargs__(self):
86
+        return (self.family, self.tag, True)
87
+
88
+    @property
89
+    def pos(self):
90
+        """Real space position of the site.
91
+
92
+        This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
93
+        """
94
+        return self.family.pos(self.tag)
95
+
96
+
97
+class SiteArray:
98
+    """An array of sites, members of a `SiteFamily`.
99
+
100
+    Parameters
101
+    ----------
102
+    family : an instance of `SiteFamily`
103
+        The 'type' of the sites.
104
+    tags : a sequence of python objects
105
+        Sequence of unique identifiers of the sites within the
106
+        site array family, typically vectors of integers.
107
+
108
+    Raises
109
+    ------
110
+    ValueError
111
+        If `tags` are not proper tags for `family`.
112
+
113
+    See Also
114
+    --------
115
+    kwant.system.Site
116
+    """
117
+
118
+    def __init__(self, family, tags):
119
+        self.family = family
120
+        try:
121
+            tags = family.normalize_tags(tags)
122
+        except (TypeError, ValueError) as e:
123
+            msg = 'Tags {0} are not allowed for site family {1}: {2}'
124
+            raise type(e)(msg.format(repr(tags), repr(family), e.args[0]))
125
+        self.tags = tags
126
+
127
+    def __repr__(self):
128
+        return 'SiteArray({0}, {1})'.format(repr(self.family), repr(self.tags))
129
+
130
+    def __str__(self):
131
+        sf = self.family
132
+        return ('<SiteArray {0} of {1}>'
133
+                .format(self.tags, sf.name if sf.name else sf))
134
+
135
+    def __len__(self):
136
+        return len(self.tags)
137
+
138
+    def __eq__(self, other):
139
+        if not isinstance(other, SiteArray):
140
+            raise NotImplementedError()
141
+        return self.family == other.family and np.all(self.tags == other.tags)
142
+
143
+    def positions(self):
144
+        """Real space position of the site.
145
+
146
+        This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
147
+        """
148
+        return self.family.positions(self.tags)
149
+
150
+
151
+@total_ordering
152
+class SiteFamily:
153
+    """Abstract base class for site families.
154
+
155
+    Site families are the 'type' of `Site` objects.  Within a family, individual
156
+    sites are uniquely identified by tags.  Valid tags must be hashable Python
157
+    objects, further details are up to the family.
158
+
159
+    Site families must be immutable and fully defined by their initial
160
+    arguments.  They must inherit from this abstract base class and call its
161
+    __init__ function providing it with two arguments: a canonical
162
+    representation and a name.  The canonical representation will be returned as
163
+    the objects representation and must uniquely identify the site family
164
+    instance.  The name is a string used to distinguish otherwise identical site
165
+    families.  It may be empty. ``norbs`` defines the number of orbitals
166
+    on sites associated with this site family; it may be `None`, in which case
167
+    the number of orbitals is not specified.
168
+
169
+
170
+    All site families must define either 'normalize_tag' or 'normalize_tags',
171
+    which brings a tag (or, in the latter case, a sequence of tags) to the
172
+    standard format for this site family.
173
+
174
+    Site families may also implement methods ``pos(tag)`` and
175
+    ``positions(tags)``, which return a vector of realspace coordinates or an
176
+    array of vectors of realspace coordinates of the site(s) belonging to this
177
+    family with the given tag(s). These methods are used in plotting routines.
178
+    ``positions(tags)`` should return an array with shape ``(N, M)`` where
179
+    ``N`` is the length of ``tags``, and ``M`` is the realspace dimension.
180
+
181
+    If the ``norbs`` of a site family are provided, and sites of this family
182
+    are used to populate a `~kwant.builder.Builder`, then the associated
183
+    Hamiltonian values must have the correct shape. That is, if a site family
184
+    has ``norbs = 2``, then any on-site terms for sites belonging to this
185
+    family should be 2x2 matrices. Similarly, any hoppings to/from sites
186
+    belonging to this family must have a matrix structure where there are two
187
+    rows/columns. This condition applies equally to Hamiltonian values that
188
+    are given by functions. If this condition is not satisfied, an error will
189
+    be raised.
190
+    """
191
+
192
+    def __init__(self, canonical_repr, name, norbs):
193
+        self.canonical_repr = canonical_repr
194
+        self.hash = hash(canonical_repr)
195
+        self.name = name
196
+        if norbs is None:
197
+            warnings.warn("Not specfying norbs is deprecated. Always specify "
198
+                          "norbs when creating site families.",
199
+                          KwantDeprecationWarning, stacklevel=3)
200
+        if norbs is not None:
201
+            if int(norbs) != norbs or norbs <= 0:
202
+                raise ValueError('The norbs parameter must be an integer > 0.')
203
+            norbs = int(norbs)
204
+        self.norbs = norbs
205
+
206
+    def __init_subclass__(cls, **kwargs):
207
+        super().__init_subclass__(**kwargs)
208
+        if (cls.normalize_tag is SiteFamily.normalize_tag
209
+            and cls.normalize_tags is SiteFamily.normalize_tags):
210
+            raise TypeError("Must redefine either 'normalize_tag' or "
211
+                            "'normalize_tags'")
212
+
213
+    def __repr__(self):
214
+        return self.canonical_repr
215
+
216
+    def __str__(self):
217
+        if self.name:
218
+            msg = '<{0} site family {1}{2}>'
219
+        else:
220
+            msg = '<unnamed {0} site family{2}>'
221
+        orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
222
+        return msg.format(self.__class__.__name__, self.name, orbs)
223
+
224
+    def __hash__(self):
225
+        return self.hash
226
+
227
+    def __eq__(self, other):
228
+        try:
229
+            return self.canonical_repr == other.canonical_repr
230
+        except AttributeError:
231
+            return False
232
+
233
+    def __ne__(self, other):
234
+        try:
235
+            return self.canonical_repr != other.canonical_repr
236
+        except AttributeError:
237
+            return True
238
+
239
+    def __lt__(self, other):
240
+        # If this raises an AttributeError, we were trying
241
+        # to compare it to something non-comparable anyway.
242
+        return self.canonical_repr < other.canonical_repr
243
+
244
+    def normalize_tag(self, tag):
245
+        """Return a normalized version of the tag.
246
+
247
+        Raises TypeError or ValueError if the tag is not acceptable.
248
+        """
249
+        tag, = self.normalize_tags([tag])
250
+        return tag
251
+
252
+    def normalize_tags(self, tags):
253
+        """Return a normalized version of the tags.
254
+
255
+        Raises TypeError or ValueError if the tags are not acceptable.
256
+        """
257
+        return np.array([self.normalize_tag(tag) for tag in tags])
258
+
259
+    def __call__(self, *tag):
260
+        """
261
+        A convenience function.
262
+
263
+        This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
264
+        """
265
+        # Catch a likely and difficult to find mistake.
266
+        if tag and isinstance(tag[0], tuple):
267
+            raise ValueError('Use site_family(1, 2) instead of '
268
+                             'site_family((1, 2))!')
269
+        return Site(self, tag)
270
+
271
+
272
+
273
+################ Systems
18 274
 
19 275
 
20 276
 class System(metaclass=abc.ABCMeta):
... ...
@@ -92,12 +348,100 @@ class System(metaclass=abc.ABCMeta):
92 348
         details = ', and '.join((', '.join(details[:-1]), details[-1]))
93 349
         return '<{} with {}>'.format(self.__class__.__name__, details)
94 350
 
351
+    hamiltonian_submatrix = _system.hamiltonian_submatrix
352
+
353
+
354
+Term = namedtuple(
355
+    "Term",
356
+    ["subgraph", "hermitian", "parameters"],
357
+)
358
+
359
+
360
+class VectorizedSystem(System, metaclass=abc.ABCMeta):
361
+    """Abstract general low-level system with support for vectorization.
362
+
363
+    Attributes
364
+    ----------
365
+    graph : kwant.graph.CGraph
366
+        The system graph.
367
+    subgraphs : sequence of tuples
368
+        Each subgraph has the form '((idx1, idx2), (offsets1, offsets2))'
369
+        where 'offsets1' and 'offsets2' index sites within the site arrays
370
+        indexed by 'idx1' and 'idx2'.
371
+    terms : sequence of tuples
372
+        Each tuple has the following structure:
373
+        (subgraph: int, hermitian: bool, parameters: List(str))
374
+        'subgraph' indexes 'subgraphs' and supplies the to/from sites of this
375
+        term. 'hermitian' is 'True' if the term needs its Hermitian
376
+        conjugate to be added when evaluating the Hamiltonian, and 'parameters'
377
+        contains a list of parameter names used when evaluating this term.
378
+    site_arrays : sequence of SiteArray
379
+        The sites of the system.
380
+    site_ranges : None or Nx3 integer array
381
+        Has 1 row per site array, plus one extra row.  Each row consists
382
+        of ``(first_site, norbs, orb_offset)``: the index of the first
383
+        site in the site array, the number of orbitals on each site in
384
+        the site array, and the offset of the first orbital of the first
385
+        site in the site array.  In addition, the final row has the form
386
+        ``(len(graph.num_nodes), 0, tot_norbs)`` where ``tot_norbs`` is the
387
+        total number of orbitals in the system.  ``None`` if any site array
388
+        in 'site_arrays' does not have 'norbs' specified. Note 'site_ranges'
389
+        is directly computable from 'site_arrays'.
390
+    parameters : frozenset of strings
391
+        The names of the parameters on which the system depends. This attribute
392
+        is provisional and may be changed in a future version of Kwant
393
+
394
+    Notes
395
+    -----
396
+    The sites of the system are indexed by integers ranging from 0 to
397
+    ``self.graph.num_nodes - 1``.
398
+
399
+    Optionally, a class derived from ``System`` can provide a method ``pos`` which
400
+    is assumed to return the real-space position of a site given its index.
401
+    """
402
+    @abc.abstractmethod
403
+    def hamiltonian_term(self, term_number, selector=slice(None),
404
+                         args=(), params=None):
405
+        """Return the Hamiltonians for hamiltonian term number k.
95 406
 
96
-# Add a C-implemented function as an unbound method to class System.
97
-System.hamiltonian_submatrix = _system.hamiltonian_submatrix
407
+        Parameters
408
+        ----------
409
+        term_number : int
410
+            The number of the term to evaluate.
411
+        selector : slice or sequence of int, default: slice(None)
412
+            The elements of the term to evaluate.
413
+        args : tuple
414
+            Positional arguments to the term. (Deprecated)
415
+        params : dict
416
+            Keyword parameters to the term
98 417
 
418
+        Returns
419
+        -------
420
+        hamiltonian : 3d complex array
421
+            Has shape ``(N, P, Q)`` where ``N`` is the number of matrix
422
+            elements in this term (or the number selected by 'selector'
423
+            if provided), ``P`` and ``Q`` are the number of orbitals in the
424
+            'to' and 'from' site arrays associated with this term.
99 425
 
100
-class FiniteSystem(System, metaclass=abc.ABCMeta):
426
+        Providing positional arguments via 'args' is deprecated,
427
+        instead, provide named parameters as a dictionary via 'params'.
428
+        """
429
+    @property
430
+    @lru_cache(1)
431
+    def site_ranges(self):
432
+        site_offsets = np.cumsum([0] + [len(arr) for arr in self.site_arrays])
433
+        norbs = [arr.family.norbs for arr in self.site_arrays] + [0]
434
+        if any(norb is None for norb in norbs):
435
+            return None
436
+        orb_offsets = np.cumsum(
437
+            [0] + [len(arr) * arr.family.norbs for arr in self.site_arrays]
438
+        )
439
+        return np.array([site_offsets, norbs, orb_offsets]).transpose()
440
+
441
+    hamiltonian_submatrix = _system.vectorized_hamiltonian_submatrix
442
+
443
+
444
+class FiniteSystemMixin(metaclass=abc.ABCMeta):
101 445
     """Abstract finite low-level system, possibly with leads.
102 446
 
103 447
     Attributes
... ...
@@ -220,7 +564,19 @@ class FiniteSystem(System, metaclass=abc.ABCMeta):
220 564
         return symmetries.validate(ham)
221 565
 
222 566
 
223
-class InfiniteSystem(System, metaclass=abc.ABCMeta):
567
+class FiniteSystem(System, FiniteSystemMixin, metaclass=abc.ABCMeta):
568
+    pass
569
+
570
+
571
+class FiniteVectorizedSystem(VectorizedSystem, FiniteSystemMixin, metaclass=abc.ABCMeta):
572
+    pass
573
+
574
+
575
+def is_finite(syst):
576
+    return isinstance(syst, (FiniteSystem, FiniteVectorizedSystem))
577
+
578
+
579
+class InfiniteSystemMixin(metaclass=abc.ABCMeta):
224 580
     """Abstract infinite low-level system.
225 581
 
226 582
     An infinite system consists of an infinite series of identical cells.
... ...
@@ -261,30 +617,10 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
261 617
     infinite system.  The other scheme has the numbers of site 0 and 1
262 618
     exchanged, as well as of site 3 and 4.
263 619
 
620
+    Sites in the fundamental domain cell must belong to a different site array
621
+    than the sites in the previous cell. In the above example this means that
622
+    sites '(0, 1, 2)' and '(3, 4)' must belong to different site arrays.
264 623
     """
265
-    @deprecate_args
266
-    def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
267
-        """Hamiltonian of a single cell of the infinite system.
268
-
269
-        Providing positional arguments via 'args' is deprecated,
270
-        instead, provide named parameters as a dictionary via 'params'.
271
-        """
272
-        cell_sites = range(self.cell_size)
273
-        return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
274
-                                          sparse=sparse, params=params)
275
-
276
-    @deprecate_args
277
-    def inter_cell_hopping(self, args=(), sparse=False, *, params=None):
278
-        """Hopping Hamiltonian between two cells of the infinite system.
279
-
280
-        Providing positional arguments via 'args' is deprecated,
281
-        instead, provide named parameters as a dictionary via 'params'.
282
-        """
283
-        cell_sites = range(self.cell_size)
284
-        interface_sites = range(self.cell_size, self.graph.num_nodes)
285
-        return self.hamiltonian_submatrix(args, cell_sites, interface_sites,
286
-                                          sparse=sparse, params=params)
287
-
288 624
     @deprecate_args
289 625
     def modes(self, energy=0, args=(), *, params=None):
290 626
         """Return mode decomposition of the lead
... ...
@@ -368,6 +704,41 @@ class InfiniteSystem(System, metaclass=abc.ABCMeta):
368 704
         return list(broken)
369 705
 
370 706
 
707
+class InfiniteSystem(System, InfiniteSystemMixin, metaclass=abc.ABCMeta):
708
+
709
+    @deprecate_args
710
+    def cell_hamiltonian(self, args=(), sparse=False, *, params=None):
711
+        """Hamiltonian of a single cell of the infinite system.
712
+
713
+        Providing positional arguments via 'args' is deprecated,
714
+        instead, provide named parameters as a dictionary via 'params'.
715
+        """
716
+        cell_sites = range(self.cell_size)
717
+        return self.hamiltonian_submatrix(args, cell_sites, cell_sites,
718
+                                          sparse=sparse, params=params)
719
+
720
+    @deprecate_args
721
+    def inter_cell_hopping(self, args=(), sparse=False, *, params=None):
722
+        """Hopping Hamiltonian between two cells of the infinite system.
723
+
724
+        Providing positional arguments via 'args' is deprecated,
725
+        instead, provide named parameters as a dictionary via 'params'.
726
+        """
727
+        cell_sites = range(self.cell_size)
728
+        interface_sites = range(self.cell_size, self.graph.num_nodes)
729
+        return self.hamiltonian_submatrix(args, cell_sites, interface_sites,
730
+                                          sparse=sparse, params=params)
731
+
732
+
733
+class InfiniteVectorizedSystem(VectorizedSystem, InfiniteSystemMixin, metaclass=abc.ABCMeta):
734
+    cell_hamiltonian = _system.vectorized_cell_hamiltonian
735
+    inter_cell_hopping = _system.vectorized_inter_cell_hopping
736
+
737
+
738
+def is_infinite(syst):
739
+    return isinstance(syst, (InfiniteSystem, InfiniteVectorizedSystem))
740
+
741
+
371 742
 class PrecalculatedLead:
372 743
     def __init__(self, modes=None, selfenergy=None):
373 744
         """A general lead defined by its self energy.
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2018 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -19,8 +19,8 @@ from pytest import raises, warns
19 19
 from numpy.testing import assert_almost_equal
20 20
 
21 21
 import kwant
22
-from kwant import builder
23
-from kwant._common import ensure_rng
22
+from kwant import builder, system
23
+from kwant._common import ensure_rng, KwantDeprecationWarning
24 24
 
25 25
 
26 26
 def test_bad_keys():
... ...
@@ -290,12 +290,28 @@ def random_hopping_integral(rng):
290 290
 
291 291
 
292 292
 def check_onsite(fsyst, sites, subset=False, check_values=True):
293
+    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
294
+
295
+    if vectorized:
296
+        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
297
+
293 298
     freq = {}
294 299
     for node in range(fsyst.graph.num_nodes):
295 300
         site = fsyst.sites[node].tag
296 301
         freq[site] = freq.get(site, 0) + 1
297 302
         if check_values and site in sites:
298
-            assert fsyst.onsites[node][0] is sites[site]
303
+            if vectorized:
304
+                term = fsyst._onsite_term_by_site_id[node]
305
+                value = fsyst._term_values[term]
306
+                if callable(value):
307
+                    assert value is sites[site]
308
+                else:
309
+                    (w, _), (off, _) = fsyst.subgraphs[fsyst.terms[term].subgraph]
310
+                    node_off = node - site_offsets[w]
311
+                    selector = np.searchsorted(off, node_off)
312
+                    assert np.allclose(value[selector], sites[site])
313
+            else:
314
+                assert fsyst.onsites[node][0] is sites[site]
299 315
     if not subset:
300 316
         # Check that all sites of `fsyst` are in `sites`.
301 317
         for site in freq.keys():
... ...
@@ -306,24 +322,50 @@ def check_onsite(fsyst, sites, subset=False, check_values=True):
306 322
 
307 323
 
308 324
 def check_hoppings(fsyst, hops):
325
+    vectorized = isinstance(fsyst, (system.FiniteVectorizedSystem, system.InfiniteVectorizedSystem))
326
+
327
+    if vectorized:
328
+        site_offsets = np.cumsum([0] + [len(s) for s in fsyst.site_arrays])
329
+
309 330
     assert fsyst.graph.num_edges == 2 * len(hops)
310 331
     for edge_id, edge in enumerate(fsyst.graph):
311
-        tail, head = edge
312
-        tail = fsyst.sites[tail].tag
313
-        head = fsyst.sites[head].tag
314
-        value = fsyst.hoppings[edge_id][0]
315
-        if value is builder.Other:
316
-            assert (head, tail) in hops
332
+        i, j = edge
333
+        tail = fsyst.sites[i].tag
334
+        head = fsyst.sites[j].tag
335
+
336
+        if vectorized:
337
+            term = fsyst._hopping_term_by_edge_id[edge_id]
338
+            if term < 0:  # Hermitian conjugate
339
+                assert (head, tail) in hops
340
+            else:
341
+                value = fsyst._term_values[term]
342
+                assert (tail, head) in hops
343
+                if callable(value):
344
+                    assert value is hops[tail, head]
345
+                else:
346
+                    dtype = np.dtype([('f0', int), ('f1', int)])
347
+                    subgraph = fsyst.terms[term].subgraph
348
+                    (to_w, from_w), hoppings = fsyst.subgraphs[subgraph]
349
+                    hop = (i - site_offsets[to_w], j - site_offsets[from_w])
350
+                    hop = np.array(hop, dtype=dtype)
351
+                    hoppings = hoppings.transpose().view(dtype).reshape(-1)
352
+                    selector = np.recarray.searchsorted(hoppings, hop)
353
+                    assert np.allclose(value[selector], hops[tail, head])
317 354
         else:
318
-            assert (tail, head) in hops
319
-            assert value is hops[tail, head]
355
+            value = fsyst.hoppings[edge_id][0]
356
+            if value is builder.Other:
357
+                assert (head, tail) in hops
358
+            else:
359
+                assert (tail, head) in hops
360
+                assert value is hops[tail, head]
320 361
 
321 362
 def check_id_by_site(fsyst):
322 363
     for i, site in enumerate(fsyst.sites):
323 364
         assert fsyst.id_by_site[site] == i
324 365
 
325 366
 
326
-def test_finalization():
367
+@pytest.mark.parametrize("vectorize", [False, True])
368
+def test_finalization(vectorize):
327 369
     """Test the finalization of finite and infinite systems.
328 370
 
329 371
     In order to exactly verify the finalization, low-level features of the
... ...
@@ -377,7 +419,7 @@ def test_finalization():
377 419
     neighbors = sorted(neighbors)
378 420
 
379 421
     # Build scattering region from blueprint and test it.
380
-    syst = builder.Builder()
422
+    syst = builder.Builder(vectorize=vectorize)
381 423
     fam = kwant.lattice.general(ta.identity(2), norbs=1)
382 424
     for site, value in sr_sites.items():
383 425
         syst[fam(*site)] = value
... ...
@@ -388,7 +430,7 @@ def test_finalization():
388 430
     check_onsite(fsyst, sr_sites)
389 431
     check_hoppings(fsyst, sr_hops)
390 432
     # check that sites are sorted
391
-    assert fsyst.sites == tuple(sorted(fam(*site) for site in sr_sites))
433
+    assert tuple(fsyst.sites) == tuple(sorted(fam(*site) for site in sr_sites))
392 434
 
393 435
     # Build lead from blueprint and test it.
394 436
     lead = builder.Builder(kwant.TranslationalSymmetry((size, 0)))
... ...
@@ -421,12 +463,12 @@ def test_finalization():
421 463
 
422 464
     # Attach lead with improper interface.
423 465
     syst.leads[-1] = builder.BuilderLead(
424
-        lead, 2 * tuple(builder.Site(fam, n) for n in neighbors))
466
+        lead, 2 * tuple(system.Site(fam, n) for n in neighbors))
425 467
     raises(ValueError, syst.finalized)
426 468
 
427 469
     # Attach lead properly.
428 470
     syst.leads[-1] = builder.BuilderLead(
429
-        lead, (builder.Site(fam, n) for n in neighbors))
471
+        lead, (system.Site(fam, n) for n in neighbors))
430 472
     fsyst = syst.finalized()
431 473
     assert len(fsyst.lead_interfaces) == 1
432 474
     assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
... ...
@@ -434,15 +476,15 @@ def test_finalization():
434 476
 
435 477
     # test that we cannot finalize a system with a badly sorted interface order
436 478
     raises(ValueError, builder.InfiniteSystem, lead,
437
-           [builder.Site(fam, n) for n in reversed(neighbors)])
479
+           [system.Site(fam, n) for n in reversed(neighbors)])
438 480
     # site ordering independent of whether interface was specified
439
-    flead_order = builder.InfiniteSystem(lead, [builder.Site(fam, n)
481
+    flead_order = builder.InfiniteSystem(lead, [system.Site(fam, n)
440 482
                                                 for n in neighbors])
441
-    assert flead.sites == flead_order.sites
483
+    assert tuple(flead.sites) == tuple(flead_order.sites)
442 484
 
443 485
 
444 486
     syst.leads[-1] = builder.BuilderLead(
445
-        lead, (builder.Site(fam, n) for n in neighbors))
487
+        lead, (system.Site(fam, n) for n in neighbors))
446 488
     fsyst = syst.finalized()
447 489
     assert len(fsyst.lead_interfaces) == 1
448 490
     assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
... ...
@@ -457,47 +499,62 @@ def test_finalization():
457 499
     raises(ValueError, lead.finalized)
458 500
 
459 501
 
460
-def test_site_ranges():
502
+def _make_system_from_sites(sites, vectorize):
503
+    syst = builder.Builder(vectorize=vectorize)
504
+    for s in sites:
505
+        norbs = s.family.norbs
506
+        if norbs:
507
+            syst[s] = np.eye(s.family.norbs)
508
+        else:
509
+            syst[s] = None
510
+    return syst.finalized()
511
+
512
+
513
+@pytest.mark.parametrize("vectorize", [False, True])
514
+def test_site_ranges(vectorize):
461 515
     lat1a = kwant.lattice.chain(norbs=1, name='a')
462 516
     lat1b = kwant.lattice.chain(norbs=1, name='b')
463 517
     lat2 = kwant.lattice.chain(norbs=2)
464
-    site_ranges = builder._site_ranges
465 518
 
466 519
     # simple case -- single site family
467 520
     for lat in (lat1a, lat2):
468 521
         sites = list(map(lat, range(10)))
469
-        ranges = site_ranges(sites)
522
+        syst = _make_system_from_sites(sites, vectorize)
523
+        ranges = syst.site_ranges
470 524
         expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
471
-        assert ranges == expected
525
+        assert np.array_equal(ranges, expected)
472 526
 
473 527
     # pair of site families
474
-    sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)),
475
-                     map(lat1a, range(4)))
476
-    expected = [(0, 1, 0), (4, 1, 4), (10, 1, 10), (14, 0, 14)]
477
-    assert expected == site_ranges(tuple(sites))
528
+    sites = it.chain(map(lat1a, range(4)), map(lat1b, range(6)))
529
+    syst = _make_system_from_sites(sites, vectorize)
530
+    expected = [(0, 1, 0), (4, 1, 4), (10, 0, 10)]
531
+    assert np.array_equal(expected, syst.site_ranges)
478 532
     sites = it.chain(map(lat2, range(4)), map(lat1a, range(6)),
479 533
                      map(lat1b, range(4)))
534
+    syst = _make_system_from_sites(sites, vectorize)
480 535
     expected = [(0, 2, 0), (4, 1, 4*2), (10, 1, 4*2+6), (14, 0, 4*2+10)]
481
-    assert expected == site_ranges(tuple(sites))
536
+    assert np.array_equal(expected, syst.site_ranges)
482 537
 
483 538
     # test with an actual builder
484 539
     for lat in (lat1a, lat2):
485 540
         sites = list(map(lat, range(10)))
486
-        syst = kwant.Builder()
541
+        syst = kwant.Builder(vectorize=vectorize)
487 542
         syst[sites] = np.eye(lat.norbs)
488 543
         ranges = syst.finalized().site_ranges
489 544
         expected = [(0, lat.norbs, 0), (10, 0, 10 * lat.norbs)]
490
-        assert ranges == expected
491
-        # poison system with a single site with no norbs defined.
492
-        # Also catch the deprecation warning.
493
-        with warnings.catch_warnings():
494
-            warnings.simplefilter("ignore")
495
-            syst[kwant.lattice.chain(norbs=None)(0)] = 1
496
-        ranges = syst.finalized().site_ranges
497
-        assert ranges == None
498
-
499
-
500
-def test_hamiltonian_evaluation():
545
+        assert np.array_equal(ranges, expected)
546
+        if not vectorize:  # vectorized systems *must* have all norbs
547
+            # poison system with a single site with no norbs defined.
548
+            # Also catch the deprecation warning.
549
+            with warnings.catch_warnings():
550
+                warnings.simplefilter("ignore")
551
+                syst[kwant.lattice.chain(norbs=None)(0)] = 1
552
+            ranges = syst.finalized().site_ranges
553
+            assert ranges is None
554
+
555
+
556
+@pytest.mark.parametrize("vectorize", [False, True])
557
+def test_hamiltonian_evaluation(vectorize):
501 558
     def f_onsite(site):
502 559
         return site.tag[0]
503 560
 
... ...
@@ -505,14 +562,26 @@ def test_hamiltonian_evaluation():
505 562
         a, b = a.tag, b.tag
506 563
         return complex(a[0] + b[0], a[1] - b[1])
507 564
 
565
+    def f_onsite_vectorized(sites):
566
+        return sites.tags[:, 0]
567
+
568
+    def f_hopping_vectorized(a, b):
569
+        a, b = a.tags, b.tags
570
+        return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])
571
+
508 572
     tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
509 573
     edges = [(0, 1), (0, 2), (0, 3), (1, 2)]
510 574
 
511
-    syst = builder.Builder()
575
+    syst = builder.Builder(vectorize=vectorize)
512 576
     fam = builder.SimpleSiteFamily(norbs=1)
513 577
     sites = [fam(*tag) for tag in tags]
514
-    syst[(fam(*tag) for tag in tags)] = f_onsite
515
-    syst[((fam(*tags[i]), fam(*tags[j])) for (i, j) in edges)] = f_hopping
578
+    hoppings = [(sites[i], sites[j]) for i, j in edges]
579
+    if vectorize:
580
+        syst[sites] = f_onsite_vectorized
581
+        syst[hoppings] = f_hopping_vectorized
582
+    else:
583
+        syst[sites] = f_onsite
584
+        syst[hoppings] = f_hopping
516 585
     fsyst = syst.finalized()
517 586
 
518 587
     assert fsyst.graph.num_nodes == len(tags)
... ...
@@ -521,12 +590,16 @@ def test_hamiltonian_evaluation():
521 590
     for i in range(len(tags)):
522 591
         site = fsyst.sites[i]
523 592
         assert site in sites
524
-        assert fsyst.hamiltonian(i, i) == syst[site](site)
593
+        with warnings.catch_warnings():
594
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
595
+            assert fsyst.hamiltonian(i, i) == f_onsite(site)
525 596
 
526 597
     for t, h in fsyst.graph:
527 598
         tsite = fsyst.sites[t]
528 599
         hsite = fsyst.sites[h]
529
-        assert fsyst.hamiltonian(t, h) == syst[tsite, hsite](tsite, hsite)
600
+        with warnings.catch_warnings():
601
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
602
+            assert fsyst.hamiltonian(t, h) == f_hopping(tsite, hsite)
530 603
 
531 604
     # test when user-function raises errors
532 605
     def onsite_raises(site):
... ...
@@ -539,13 +612,18 @@ def test_hamiltonian_evaluation():
539 612
         a, b = hop
540 613
         # exceptions are converted to kwant.UserCodeError and we add our message
541 614
         with raises(kwant.UserCodeError) as exc_info:
542
-            fsyst.hamiltonian(a, a)
615
+            with warnings.catch_warnings():
616
+                warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
617
+                fsyst.hamiltonian(a, a)
543 618
         msg = 'Error occurred in user-supplied value function "onsite_raises"'
544 619
         assert msg in str(exc_info.value)
545 620
 
546 621
         for hop in [(a, b), (b, a)]:
547 622
             with raises(kwant.UserCodeError) as exc_info:
548
-                fsyst.hamiltonian(*hop)
623
+                with warnings.catch_warnings():
624
+                    # Ignore deprecation warnings for 'hamiltonian'
625
+                    warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
626
+                    fsyst.hamiltonian(*hop)
549 627
             msg = ('Error occurred in user-supplied '
550 628
                    'value function "hopping_raises"')
551 629
             assert msg in str(exc_info.value)
... ...
@@ -555,7 +633,7 @@ def test_hamiltonian_evaluation():
555 633
     syst[new_hop[0]] = onsite_raises
556 634
     syst[new_hop] = hopping_raises
557 635
     fsyst = syst.finalized()
558
-    hop = tuple(map(fsyst.sites.index, new_hop))
636
+    hop = tuple(map(fsyst.id_by_site.__getitem__, new_hop))
559 637
     test_raising(fsyst, hop)
560 638
 
561 639
     # test with infinite system
... ...
@@ -563,10 +641,118 @@ def test_hamiltonian_evaluation():
563 641
     for k, v in it.chain(syst.site_value_pairs(), syst.hopping_value_pairs()):
564 642
         inf_syst[k] = v
565 643
     inf_fsyst = inf_syst.finalized()
566
-    hop = tuple(map(inf_fsyst.sites.index, new_hop))
644
+    hop = tuple(map(inf_fsyst.id_by_site.__getitem__, new_hop))
567 645
     test_raising(inf_fsyst, hop)
568 646
 
569 647
 
648
+def test_vectorized_hamiltonian_evaluation():
649
+
650
+    def onsite(site):
651
+        return site.tag[0]
652
+
653
+    def vectorized_onsite(sites):
654
+        return sites.tags[:, 0]
655
+
656
+    def hopping(to_site, from_site):
657
+        a, b = to_site.tag, from_site.tag
658
+        return a[0] + b[0] + 1j * (a[1] - b[1])
659
+
660
+    def vectorized_hopping(to_sites, from_sites):
661
+        a, b = to_sites.tags, from_sites.tags
662
+        return a[:, 0] + b[:, 0] + 1j * (a[:, 1] - b[:, 1])
663
+
664
+    tags = [(0, 0), (1, 1), (2, 2), (3, 3)]
665
+    edges = [(0, 1), (0, 2), (0, 3), (1, 2)]
666
+
667
+    fam = builder.SimpleSiteFamily(norbs=1)
668
+    sites = [fam(*tag) for tag in tags]
669
+    hops = [(fam(*tags[i]), fam(*tags[j])) for (i, j) in edges]
670
+
671
+    syst_simple = builder.Builder(vectorize=False)
672
+    syst_simple[sites] = onsite
673
+    syst_simple[hops] = hopping
674
+    fsyst_simple = syst_simple.finalized()
675
+
676
+    syst_vectorized = builder.Builder(vectorize=True)
677
+    syst_vectorized[sites] = vectorized_onsite
678
+    syst_vectorized[hops] = vectorized_hopping
679
+    fsyst_vectorized = syst_vectorized.finalized()
680
+
681
+    assert fsyst_vectorized.graph.num_nodes == len(tags)
682
+    assert fsyst_vectorized.graph.num_edges == 2 * len(edges)
683
+    assert len(fsyst_vectorized.site_arrays) == 1
684
+    assert fsyst_vectorized.site_arrays[0] == system.SiteArray(fam, tags)
685
+
686
+    assert np.allclose(
687
+        fsyst_simple.hamiltonian_submatrix(),
688
+        fsyst_vectorized.hamiltonian_submatrix(),
689
+    )
690
+
691
+    for i in range(len(tags)):
692
+        site = fsyst_vectorized.sites[i]
693
+        assert site in sites
694
+        with warnings.catch_warnings():
695
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
696
+            assert (
697
+                fsyst_vectorized.hamiltonian(i, i)
698
+                == fsyst_simple.hamiltonian(i, i))
699
+
700
+    for t, h in fsyst_vectorized.graph:
701
+        with warnings.catch_warnings():
702
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
703
+            assert (
704
+                fsyst_vectorized.hamiltonian(t, h)
705
+                == fsyst_simple.hamiltonian(t, h))
706
+
707
+    # Test infinite system, including hoppings that go both ways into
708
+    # the next cell
709
+    lat = kwant.lattice.square(norbs=1)
710
+
711
+    syst_vectorized = builder.Builder(kwant.TranslationalSymmetry((-1, 0)),
712
+                                      vectorize=True)
713
+    syst_vectorized[lat(0, 0)] = 4
714
+    syst_vectorized[lat(0, 1)] = 5
715
+    syst_vectorized[lat(0, 2)] = vectorized_onsite
716
+    syst_vectorized[(lat(1, 0), lat(0, 0))] = 1j
717
+    syst_vectorized[(lat(2, 1), lat(1, 1))] = vectorized_hopping
718
+    fsyst_vectorized = syst_vectorized.finalized()
719
+
720
+    syst_simple = builder.Builder(kwant.TranslationalSymmetry((-1, 0)),
721
+                                      vectorize=False)
722
+    syst_simple[lat(0, 0)] = 4
723
+    syst_simple[lat(0, 1)] = 5
724
+    syst_simple[lat(0, 2)] = onsite
725
+    syst_simple[(lat(1, 0), lat(0, 0))] = 1j
726
+    syst_simple[(lat(2, 1), lat(1, 1))] = hopping
727
+    fsyst_simple = syst_simple.finalized()
728
+
729
+    assert np.allclose(
730
+        fsyst_vectorized.hamiltonian_submatrix(),
731
+        fsyst_simple.hamiltonian_submatrix(),
732
+    )
733
+    assert np.allclose(
734
+        fsyst_vectorized.cell_hamiltonian(),
735
+        fsyst_simple.cell_hamiltonian(),
736
+    )
737
+    assert np.allclose(
738
+        fsyst_vectorized.inter_cell_hopping(),
739
+        fsyst_simple.inter_cell_hopping(),
740
+    )
741
+
742
+
743
+def test_vectorized_requires_norbs():
744
+
745
+    # Catch deprecation warning for lack of norbs
746
+    with warnings.catch_warnings():
747
+        warnings.simplefilter("ignore")
748
+        fam = builder.SimpleSiteFamily()
749
+
750
+    syst = builder.Builder(vectorize=True)
751
+    syst[fam(0, 0)] = 1
752
+
753
+    raises(ValueError, syst.finalized)
754
+
755
+
570 756
 def test_dangling():
571 757
     def make_system():
572 758
         #        1
... ...
@@ -1204,15 +1390,22 @@ def test_discrete_symmetries():
1204 1390
 # We need to keep testing 'args', but we don't want to see
1205 1391
 # all the deprecation warnings in the test logs
1206 1392
 @pytest.mark.filterwarnings("ignore:.*'args' parameter")
1207
-def test_argument_passing():
1393
+@pytest.mark.parametrize("vectorize", [False, True])
1394
+def test_argument_passing(vectorize):
1208 1395
     chain = kwant.lattice.chain(norbs=1)
1209 1396
 
1210 1397
     # Test for passing parameters to hamiltonian matrix elements
1211 1398
     def onsite(site, p1, p2):
1212
-        return p1 + p2
1399
+        if vectorize:
1400
+            return (p1 + p2) * np.ones(len(site))
1401
+        else:
1402
+            return p1 + p2
1213 1403
 
1214 1404
     def hopping(site1, site2, p1, p2):
1215
-        return p1 - p2
1405
+        if vectorize:
1406
+            return (p1 - p2) * np.ones(len(site1))
1407
+        else:
1408
+            return p1 - p2
1216 1409
 
1217 1410
     def gen_fill_syst(onsite, hopping, syst):
1218 1411
         syst[(chain(i) for i in range(3))] = onsite
... ...
@@ -1221,8 +1414,9 @@ def test_argument_passing():
1221 1414
 
1222 1415
     fill_syst = ft.partial(gen_fill_syst, onsite, hopping)
1223 1416
 
1224
-    syst = fill_syst(kwant.Builder())
1225
-    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
1417
+    syst = fill_syst(kwant.Builder(vectorize=vectorize))
1418
+    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)),
1419
+                                       vectorize=vectorize))
1226 1420
 
1227 1421
     tests = (
1228 1422
         syst.hamiltonian_submatrix,
... ...
@@ -1238,28 +1432,34 @@ def test_argument_passing():
1238 1432
 
1239 1433
     # test that mixing 'args' and 'params' raises TypeError
1240 1434
     with raises(TypeError):
1241
-        syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
1435
+        with warnings.catch_warnings():
1436
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
1437
+            syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
1438
+
1242 1439
     with raises(TypeError):
1243
-        inf_syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
1440
+        with warnings.catch_warnings():
1441
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
1442
+            inf_syst.hamiltonian(0, 0, *(2, 1), params=dict(p1=2, p2=1))
1244 1443
 
1245 1444
     # Missing parameters raises TypeError
1246 1445
     with raises(TypeError):
1247
-        syst.hamiltonian(0, 0, params=dict(p1=2))
1446
+        with warnings.catch_warnings():
1447
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
1448
+            syst.hamiltonian(0, 0, params=dict(p1=2))
1449
+
1248 1450
     with raises(TypeError):
1249
-        syst.hamiltonian_submatrix(params=dict(p1=2))
1451
+        with warnings.catch_warnings():
1452
+            warnings.filterwarnings("ignore", category=KwantDeprecationWarning)
1453
+            syst.hamiltonian_submatrix(params=dict(p1=2))
1250 1454
 
1251 1455
     # test that passing parameters without default values works, and that
1252 1456
     # passing parameters with default values fails
1253
-    def onsite(site, p1, p2):
1254
-        return p1 + p2
1255
-
1256
-    def hopping(site, site2, p1, p2):
1257
-        return p1 - p2
1258 1457
 
1259 1458
     fill_syst = ft.partial(gen_fill_syst, onsite, hopping)
1260 1459
 
1261
-    syst = fill_syst(kwant.Builder())
1262
-    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,))))
1460
+    syst = fill_syst(kwant.Builder(vectorize=vectorize))
1461
+    inf_syst = fill_syst(kwant.Builder(kwant.TranslationalSymmetry((-3,)),
1462
+                                       vectorize=vectorize))
1263 1463
 
1264 1464
     tests = (
1265 1465
         syst.hamiltonian_submatrix,
... ...
@@ -1275,12 +1475,18 @@ def test_argument_passing():
1275 1475
 
1276 1476
     # Some common, some different args for value functions
1277 1477
     def onsite2(site, a, b):
1278
-        return site.pos + a + b
1478
+        if vectorize:
1479
+            return site.positions()[:, 0] + a + b
1480
+        else:
1481
+            return site.pos[0] + a + b
1279 1482
 
1280 1483
     def hopping2(site1, site2, a, c, b):
1281
-        return a + b + c
1484
+        if vectorize:
1485
+            return (a + b + c) * np.ones(len(site1))
1486
+        else:
1487
+            return a + b + c
1282 1488
 
1283
-    syst = kwant.Builder()
1489
+    syst = kwant.Builder(vectorize=vectorize)
1284 1490
     syst[(chain(i) for i in range(3))] = onsite2
1285 1491
     syst[((chain(i), chain(i + 1)) for i in range(2))] = hopping2
1286 1492
     fsyst = syst.finalized()
... ...
@@ -1391,6 +1597,7 @@ def test_subs():
1391 1597
     lead = lead.finalized()
1392 1598
     assert lead.parameters == {'lead_a', 'lead_b', 'lead_c'}
1393 1599
 
1600
+
1394 1601
 def test_attach_stores_padding():
1395 1602
     lat = kwant.lattice.chain(norbs=1)
1396 1603
     syst = kwant.Builder()
... ...
@@ -445,7 +445,7 @@ def test_arg_passing(A):
445 445
     lat1 = kwant.lattice.chain(norbs=1)
446 446
 
447 447
     syst = kwant.Builder()
448
-    syst[lat1(0)] = syst[lat1(1)] = lambda s0, a, b: s0.pos + a + b
448
+    syst[lat1(0)] = syst[lat1(1)] = lambda s0, a, b: s0.pos[0] + a + b
449 449
     syst[lat1.neighbors()] = lambda s0, s1, a, b: a - b
450 450
     fsyst = syst.finalized()
451 451
 
... ...
@@ -1,4 +1,4 @@
1
-# Copyright 2011-2016 Kwant authors.
1
+# Copyright 2011-2019 Kwant authors.
2 2
 #
3 3
 # This file is part of Kwant.  It is subject to the license terms in the file
4 4
 # LICENSE.rst found in the top-level directory of this distribution and at
... ...
@@ -8,6 +8,7 @@
8 8
 
9 9
 import pickle
10 10
 import copy
11
+import pytest
11 12
 from pytest import raises
12 13
 import numpy as np
13 14
 from scipy import sparse
... ...
@@ -15,9 +16,11 @@ import kwant
15 16
 from kwant._common import ensure_rng
16 17
 
17 18
 
18
-def test_hamiltonian_submatrix():
19
-    syst = kwant.Builder()
19
+@pytest.mark.parametrize("vectorize", [False, True])
20
+def test_hamiltonian_submatrix(vectorize):
21
+    syst = kwant.Builder(vectorize=vectorize)
20 22
     chain = kwant.lattice.chain(norbs=1)
23
+    chain2 = kwant.lattice.chain(norbs=2)
21 24
     for i in range(3):
22 25
         syst[chain(i)] = 0.5 * i
23 26
     for i in range(2):
... ...
@@ -27,7 +30,11 @@ def test_hamiltonian_submatrix():
27 30
     mat = syst2.hamiltonian_submatrix()
28 31
     assert mat.shape == (3, 3)
29 32
     # Sorting is required due to unknown compression order of builder.
30
-    perm = np.argsort([os[0] for os in syst2.onsites])
33
+    if vectorize:
34
+        _, (site_offsets, _) = syst2.subgraphs[0]
35
+    else:
36
+        site_offsets = [os[0] for os in syst2.onsites]
37
+    perm = np.argsort(site_offsets)
31 38
     mat_should_be = np.array([[0, 1j, 0], [-1j, 0.5, 2j], [0, -2j, 1]])
32 39
 
33 40
     mat = mat[perm, :]
... ...
@@ -41,19 +48,12 @@ def test_hamiltonian_submatrix():
41 48
     mat = mat[:, perm]
42 49
     np.testing.assert_array_equal(mat, mat_should_be)
43 50
 
44
-    mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]])
45
-    np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
46
-
47
-    mat = syst2.hamiltonian_submatrix((), perm[[0, 1]], perm[[2]], sparse=True)
48
-    mat = mat.toarray()
49
-    np.testing.assert_array_equal(mat, mat_should_be[:2, 2:3])
50
-
51 51
     # Test for correct treatment of matrix input.
52
-    syst = kwant.Builder()
53
-    syst[chain(0)] = np.array([[0, 1j], [-1j, 0]])
52
+    syst = kwant.Builder(vectorize=vectorize)
53
+    syst[chain2(0)] = np.array([[0, 1j], [-1j, 0]])
54 54
     syst[chain(1)] = np.array([[1]])
55 55
     syst[chain(2)] = np.array([[2]])
56
-    syst[chain(1), chain(0)] = np.array([[1, 2j]])
56
+    syst[chain(1), chain2(0)] = np.array([[1, 2j]])
57 57
     syst[chain(2), chain(1)] = np.array([[3j]])
58 58
     syst2 = syst.finalized()
59 59
     mat_dense = syst2.hamiltonian_submatrix()
... ...
@@ -62,9 +62,10 @@ def test_hamiltonian_submatrix():
62 62
 
63 63
     # Test precalculation of modes.
64 64
     rng = ensure_rng(5)
65
-    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)))
66
-    lead[chain(0)] = np.zeros((2, 2))
67
-    lead[chain(0), chain(1)] = rng.randn(2, 2)
65
+    lead = kwant.Builder(kwant.TranslationalSymmetry((-1,)),
66
+                         vectorize=vectorize)
67
+    lead[chain2(0)] = np.zeros((2, 2))
68
+    lead[chain2(0), chain2(1)] = rng.randn(2, 2)
68 69
     syst.attach_lead(lead)
69 70
     syst2 = syst.finalized()
70 71
     smatrix = kwant.smatrix(syst2, .1).data
... ...
@@ -74,19 +75,26 @@ def test_hamiltonian_submatrix():
74 75
     raises(ValueError, kwant.solvers.default.greens_function, syst3, 0.2)
75 76
 
76 77
     # Test for shape errors.
77
-    syst[chain(0), chain(2)] = np.array([[1, 2]])
78
+    syst[chain2(0), chain(2)] = np.array([[1, 2]])
78 79
     syst2 = syst.finalized()
79 80
     raises(ValueError, syst2.hamiltonian_submatrix)
80 81
     raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
81
-    syst[chain(0), chain(2)] = 1
82
+    syst[chain2(0), chain(2)] = 1
82 83
     syst2 = syst.finalized()
83 84
     raises(ValueError, syst2.hamiltonian_submatrix)
84 85
     raises(ValueError, syst2.hamiltonian_submatrix, sparse=True)
86
+    if vectorize:  # non-vectorized systems don't check this at finalization
87
+        # Add another hopping of the same type but with a different
88
+        # (and still incompatible) shape.
89
+        syst[chain2(0), chain(1)] = np.array([[1, 2]])
90
+        raises(ValueError, syst.finalized)
85 91
 
86 92
 
87
-def test_pickling():
88
-    syst = kwant.Builder()
89
-    lead = kwant.Builder(symmetry=kwant.TranslationalSymmetry([1.]))
93
+@pytest.mark.parametrize("vectorize", [False, True])
94
+def test_pickling(vectorize):
95
+    syst = kwant.Builder(vectorize=vectorize)
96
+    lead = kwant.Builder(symmetry=kwant.TranslationalSymmetry([1.]),
97
+                         vectorize=vectorize)
90 98
     lat = kwant.lattice.chain(norbs=1)
91 99
     syst[lat(0)] = syst[lat(1)] = 0
92 100
     syst[lat(0), lat(1)] = 1
... ...
@@ -183,6 +183,9 @@ def wraparound(builder, keep=None, *, coordinate_names='xyz'):
183 183
         f.__signature__ = inspect.Signature(params.values())
184 184
         return f
185 185
 
186
+    if builder.vectorize:
187
+        raise TypeError("'wraparound' does not work with vectorized Builders.")
188
+
186 189
     try:
187 190
         momenta = ['k_{}'.format(coordinate_names[i])
188 191
                    for i in range(len(builder.symmetry.periods))]