diff --git a/kwant/builder.py b/kwant/builder.py
index 1215e1991a0118ed61dac38b46794101d1e81426..bd86c2531e0ac160e0801082083aae96f7690361 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -454,7 +454,7 @@ class SelfEnergyLead(Lead):
     interface : sequence of `Site` instances
     """
     def __init__(self, selfenergy_func, interface):
-        self._selfenergy_func = selfenergy_func
+        self.selfenergy_func = selfenergy_func
         self.interface = tuple(interface)
 
     def finalized(self):
@@ -462,7 +462,7 @@ class SelfEnergyLead(Lead):
         return self
 
     def selfenergy(self, energy, args=()):
-        return self._selfenergy_func(energy, args)
+        return self.selfenergy_func(energy, args)
 
 
 class ModesLead(Lead):
diff --git a/kwant/lattice.py b/kwant/lattice.py
index 470ff6885bd23bac57e5dee771f8ca23da4f2ac3..f415cec2a599ff177a1c198588f8fc29d4ebc4b3 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -98,8 +98,8 @@ class Polyatomic(object):
         # Sequence of primitive vectors of the lattice.
         self._prim_vecs = prim_vecs
         # Precalculation of auxiliary arrays for real space calculations.
-        self._reduced_vecs, self._transf = lll.lll(prim_vecs)
-        self._voronoi = ta.dot(lll.voronoi(self._reduced_vecs), self._transf)
+        self.reduced_vecs, self.transf = lll.lll(prim_vecs)
+        self.voronoi = ta.dot(lll.voronoi(self.reduced_vecs), self.transf)
 
     def shape(self, function, start):
         """Return a key for all the lattice sites inside a given shape.
@@ -158,7 +158,7 @@ class Polyatomic(object):
                 raise ValueError('Dimensionality of start position does not '
                                  'match the space dimensionality.')
             lats = self.sublattices
-            deltas = list(self._voronoi)
+            deltas = list(self.voronoi)
 
             #### Flood-fill ####
             sites = []
@@ -421,11 +421,11 @@ class Monatomic(builder.SiteFamily, Polyatomic):
         self.offset = offset
 
         # Precalculation of auxiliary arrays for real space calculations.
-        self._reduced_vecs, self._transf = lll.lll(prim_vecs)
-        self._voronoi = ta.dot(lll.voronoi(self._reduced_vecs), self._transf)
+        self.reduced_vecs, self.transf = lll.lll(prim_vecs)
+        self.voronoi = ta.dot(lll.voronoi(self.reduced_vecs), self.transf)
 
         self.dim = dim
-        self._lattice_dim = len(prim_vecs)
+        self.lattice_dim = len(prim_vecs)
 
         if name != '':
             msg = "Monatomic lattice {0}, vectors {1}, origin {2}"
@@ -442,7 +442,7 @@ class Monatomic(builder.SiteFamily, Polyatomic):
 
     def normalize_tag(self, tag):
         tag = ta.array(tag, int)
-        if len(tag) != self._lattice_dim:
+        if len(tag) != self.lattice_dim:
             raise ValueError("Dimensionality mismatch.")
         return tag
 
@@ -455,8 +455,8 @@ class Monatomic(builder.SiteFamily, Polyatomic):
             An array with sites coordinates.
         """
         # TODO (Anton): transform to tinyarrays, once ta indexing is better.
-        return np.dot(lll.cvp(pos - self.offset, self._reduced_vecs, n),
-                      self._transf.T)
+        return np.dot(lll.cvp(pos - self.offset, self.reduced_vecs, n),
+                      self.transf.T)
 
     def closest(self, pos):
         """
@@ -497,7 +497,6 @@ class TranslationalSymmetry(builder.Symmetry):
     """
     def __init__(self, *periods):
         self.periods = ta.array(periods)
-        self._periods = self.periods
         if self.periods.ndim != 2:
             # TODO: remove the second part of the following message once
             # everybody got used to it.
@@ -539,18 +538,18 @@ class TranslationalSymmetry(builder.Symmetry):
         ValueError
             If lattice `fam` is incompatible with given periods.
         """
-        dim = self._periods.shape[1]
+        dim = self.periods.shape[1]
         if fam in self.site_family_data:
             raise KeyError('Family already processed, delete it from '
                            'site_family_data first.')
         inv = np.linalg.pinv(fam.prim_vecs)
-        bravais_periods = np.dot(self._periods, inv)
+        bravais_periods = np.dot(self.periods, inv)
         # Absolute tolerance is correct in the following since we want an error
         # relative to the closest integer.
         if not np.allclose(bravais_periods, np.round(bravais_periods),
                            rtol=0, atol=1e-8) or \
            not np.allclose([fam.vec(i) for i in bravais_periods],
-                           self._periods):
+                           self.periods):
             msg = 'Site family {0} does not have commensurate periods with ' +\
                   'symmetry {1}.'
             raise ValueError(msg.format(fam, self))
@@ -655,7 +654,7 @@ class TranslationalSymmetry(builder.Symmetry):
         The resulting symmetry has all the period vectors opposite to the
         original and an identical fundamental domain.
         """
-        result = TranslationalSymmetry(*self._periods)
+        result = TranslationalSymmetry(*self.periods)
         result.site_family_data = self.site_family_data
         result.is_reversed = not self.is_reversed
         result.periods = -self.periods
diff --git a/kwant/linalg/mumps.py b/kwant/linalg/mumps.py
index 79eaa949169cd01afa1ab51f45e8658440fa917a..170a2523df42d534c9fe805ae52a51775db07e28 100644
--- a/kwant/linalg/mumps.py
+++ b/kwant/linalg/mumps.py
@@ -23,8 +23,6 @@ orderings = { 'amd' : 0, 'amf' : 2, 'scotch' : 3, 'pord' : 4, 'metis' : 5,
 ordering_name = [ 'amd', 'user-defined', 'amf',
                   'scotch', 'pord', 'metis', 'qamd']
 
-_possible_orderings = None
-
 
 def possible_orderings():
     """Return the ordering options that are available in the current
@@ -41,13 +39,12 @@ def possible_orderings():
        A list of installed orderings that can be used in the `ordering` option
        of MUMPS.
     """
-    global _possible_orderings
 
-    if not _possible_orderings:
+    if not possible_orderings.cached:
         # Try all orderings on a small test matrix, and check which one was
         # actually used.
 
-        _possible_orderings = ['auto']
+        possible_orderings.cached = ['auto']
         for ordering in [0, 2, 3, 4, 5, 6]:
             data = np.asfortranarray([1, 1], dtype=np.complex128)
             row = np.asfortranarray([1, 2], dtype=_mumps.int_dtype)
@@ -60,9 +57,11 @@ def possible_orderings():
             instance.call()
 
             if instance.infog[7] == ordering:
-                _possible_orderings.append(ordering_name[ordering])
+                possible_orderings.cached.append(ordering_name[ordering])
+
+    return possible_orderings.cached
 
-    return _possible_orderings
+possible_orderings.cached = None
 
 
 error_messages = {
diff --git a/kwant/linalg/tests/test_mumps.py b/kwant/linalg/tests/test_mumps.py
index a316be12d0af43c3b6f8545424dceebc682f0ad5..597b24d35351038c38f7646a006e9d3cc4351de9 100644
--- a/kwant/linalg/tests/test_mumps.py
+++ b/kwant/linalg/tests/test_mumps.py
@@ -8,9 +8,9 @@
 
 try:
     from kwant.linalg.mumps import MUMPSContext, schur_complement
-    _no_mumps = False
+    no_mumps = False
 except ImportError:
-    _no_mumps = True
+    no_mumps = True
 
 from kwant.lattice import honeycomb
 from kwant.builder import Builder, HoppingKind
@@ -19,7 +19,7 @@ import numpy as np
 import scipy.sparse as sp
 from _test_utils import _Random, assert_array_almost_equal
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_lu_with_dense():
     def _test_lu_with_dense(dtype):
         rand = _Random()
@@ -48,7 +48,7 @@ def test_lu_with_dense():
     _test_lu_with_dense(np.complex128)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_schur_complement_with_dense():
     def _test_schur_complement_with_dense(dtype):
         rand = _Random()
@@ -60,7 +60,7 @@ def test_schur_complement_with_dense():
     _test_schur_complement_with_dense(np.complex128)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_error_minus_9(r=10):
     """Test if MUMPSError -9 is properly caught by increasing memory"""
 
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 1a5bf39a19f4976a6b25639577dd013e6fdd241e..8bc192a15e274f510c7ebf63398828a504326753 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -68,7 +68,7 @@ if mpl_enabled:
             self.reflen = reflen
 
         def set_linewidths(self, linewidths):
-            self._linewidths_orig = nparray_if_array(linewidths)
+            self.linewidths_orig = nparray_if_array(linewidths)
 
         def draw(self, renderer):
             if self.reflen is not None:
@@ -79,7 +79,7 @@ if mpl_enabled:
             else:
                 factor = 1
 
-            super(LineCollection, self).set_linewidths(self._linewidths_orig *
+            super(LineCollection, self).set_linewidths(self.linewidths_orig *
                                                        factor)
             return super(LineCollection, self).draw(renderer)
 
@@ -89,13 +89,13 @@ if mpl_enabled:
             super(PathCollection, self).__init__(paths, sizes=sizes, **kwargs)
 
             self.reflen = reflen
-            self._linewidths_orig = nparray_if_array(self.get_linewidths())
+            self.linewidths_orig = nparray_if_array(self.get_linewidths())
 
-            self._transforms = [matplotlib.transforms.Affine2D().scale(x) for x
+            self.transforms = [matplotlib.transforms.Affine2D().scale(x) for x
                                 in sizes]
 
         def get_transforms(self):
-            return self._transforms
+            return self.transforms
 
         def get_transform(self):
             Affine2D = matplotlib.transforms.Affine2D
@@ -112,7 +112,7 @@ if mpl_enabled:
                 # Note: only works for aspect ratio 1!
                 factor = (self.axes.transData.frozen().to_values()[0] /
                           self.figure.dpi * 72.0 * self.reflen)
-                self.set_linewidths(self._linewidths_orig * factor)
+                self.set_linewidths(self.linewidths_orig * factor)
 
             return collections.Collection.draw(self, renderer)
 
@@ -156,10 +156,10 @@ if mpl_enabled:
             def __init__(self, segments, reflen=None, zorder=0, **kwargs):
                 super(Line3DCollection, self).__init__(segments, **kwargs)
                 self.reflen = reflen
-                self._zorder3d = zorder
+                self.zorder3d = zorder
 
             def set_linewidths(self, linewidths):
-                self._linewidths_orig = nparray_if_array(linewidths)
+                self.linewidths_orig = nparray_if_array(linewidths)
 
             def do_3d_projection(self, renderer):
                 super(Line3DCollection, self).do_3d_projection(renderer)
@@ -168,7 +168,7 @@ if mpl_enabled:
                 # "-" due to the different logic in the 3d plotting, we still
                 # want larger zorder values to be plotted on top of smaller
                 # ones.
-                return -self._zorder3d
+                return -self.zorder3d
 
             def draw(self, renderer):
                 if self.reflen:
@@ -185,7 +185,7 @@ if mpl_enabled:
                     factor = 1
 
                 super(Line3DCollection, self).set_linewidths(
-                                                self._linewidths_orig * factor)
+                                                self.linewidths_orig * factor)
                 super(Line3DCollection, self).draw(renderer)
 
 
@@ -203,38 +203,38 @@ if mpl_enabled:
                     self.set_3d_properties(zs=offsets[:, 2], zdir="z")
 
                 self.reflen = reflen
-                self._zorder3d = zorder
+                self.zorder3d = zorder
 
-                self._paths_orig = np.array(paths, dtype='object')
-                self._linewidths_orig = nparray_if_array(self.get_linewidths())
-                self._linewidths_orig2 = self._linewidths_orig
-                self._array_orig = nparray_if_array(self.get_array())
-                self._facecolors_orig = nparray_if_array(self.get_facecolors())
-                self._edgecolors_orig = nparray_if_array(self.get_edgecolors())
+                self.paths_orig = np.array(paths, dtype='object')
+                self.linewidths_orig = nparray_if_array(self.get_linewidths())
+                self.linewidths_orig2 = self.linewidths_orig
+                self.array_orig = nparray_if_array(self.get_array())
+                self.facecolors_orig = nparray_if_array(self.get_facecolors())
+                self.edgecolors_orig = nparray_if_array(self.get_edgecolors())
 
                 Affine2D = matplotlib.transforms.Affine2D
-                self._orig_transforms = np.array([Affine2D().scale(x) for x in
+                self.orig_transforms = np.array([Affine2D().scale(x) for x in
                                                   sizes], dtype='object')
-                self._transforms = self._orig_transforms
+                self.transforms = self.orig_transforms
 
             def set_array(self, array):
-                self._array_orig = nparray_if_array(array)
+                self.array_orig = nparray_if_array(array)
                 super(Path3DCollection, self).set_array(array)
 
             def set_color(self, colors):
-                self._facecolors_orig = nparray_if_array(colors)
-                self._edgecolors_orig = self._facecolors_orig
+                self.facecolors_orig = nparray_if_array(colors)
+                self.edgecolors_orig = self.facecolors_orig
                 super(Path3DCollection, self).set_color(colors)
 
             def set_edgecolors(self, colors):
                 colors = matplotlib.colors.colorConverter.to_rgba_array(colors)
-                self._edgecolors_orig = nparray_if_array(colors)
+                self.edgecolors_orig = nparray_if_array(colors)
                 super(Path3DCollection, self).set_edgecolors(colors)
 
             def get_transforms(self):
                 # this is exact only for an isometric projection, for the
                 # perspective projection used in mplot3d it's an approximation
-                return self._transforms
+                return self.transforms
 
             def get_transform(self):
                 Affine2D = matplotlib.transforms.Affine2D
@@ -253,7 +253,7 @@ if mpl_enabled:
 
                 # numpy complains about zero-length index arrays
                 if len(xs) == 0:
-                    return -self._zorder3d
+                    return -self.zorder3d
 
                 proj = mplot3d.proj3d.proj_transform_clip
                 vs = np.array(proj(xs, ys, zs, renderer.M)[:3])
@@ -263,41 +263,41 @@ if mpl_enabled:
 
                     self.set_offsets(vs[:2, indx].T)
 
-                    if len(self._paths_orig) > 1:
-                        paths = np.resize(self._paths_orig, (vs.shape[1],))
+                    if len(self.paths_orig) > 1:
+                        paths = np.resize(self.paths_orig, (vs.shape[1],))
                         self.set_paths(paths[indx])
 
-                    if len(self._orig_transforms) > 1:
-                        self._transforms = np.resize(self._orig_transforms,
+                    if len(self.orig_transforms) > 1:
+                        self.transforms = np.resize(self.orig_transforms,
                                                      (vs.shape[1],))
-                        self._transforms = self._transforms[indx]
+                        self.transforms = self.transforms[indx]
 
-                    lw_orig = self._linewidths_orig
+                    lw_orig = self.linewidths_orig
                     if (isinstance(lw_orig, np.ndarray) and len(lw_orig) > 1):
-                        self._linewidths_orig2 = np.resize(lw_orig,
+                        self.linewidths_orig2 = np.resize(lw_orig,
                                                            (vs.shape[1],))[indx]
 
                     # Note: here array, facecolors and edgecolors are
                     #       guaranteed to be 2d numpy arrays or None.  (And
                     #       array is the same length as the coordinates)
 
-                    if self._array_orig is not None:
+                    if self.array_orig is not None:
                         super(Path3DCollection,
-                              self).set_array(self._array_orig[indx])
+                              self).set_array(self.array_orig[indx])
 
-                    if (self._facecolors_orig is not None and
-                        self._facecolors_orig.shape[0] > 1):
-                        shape = list(self._facecolors_orig.shape)
+                    if (self.facecolors_orig is not None and
+                        self.facecolors_orig.shape[0] > 1):
+                        shape = list(self.facecolors_orig.shape)
                         shape[0] = vs.shape[1]
                         super(Path3DCollection, self).set_facecolors(
-                            np.resize(self._facecolors_orig, shape)[indx])
+                            np.resize(self.facecolors_orig, shape)[indx])
 
-                    if (self._edgecolors_orig is not None and
-                        self._edgecolors_orig.shape[0] > 1):
-                        shape = list(self._edgecolors_orig.shape)
+                    if (self.edgecolors_orig is not None and
+                        self.edgecolors_orig.shape[0] > 1):
+                        shape = list(self.edgecolors_orig.shape)
                         shape[0] = vs.shape[1]
                         super(Path3DCollection, self).set_edgecolors(
-                                                np.resize(self._edgecolors_orig,
+                                                np.resize(self.edgecolors_orig,
                                                           shape)[indx])
                 else:
                     self.set_offsets(vs[:2].T)
@@ -316,7 +316,7 @@ if mpl_enabled:
                 proj = mplot3d.proj3d.proj_transform_clip
                 cz = proj(*(list(np.dot(corners, bbox)) + [renderer.M]))[2]
 
-                return -self._zorder3d + vs[2].mean() / cz.ptp()
+                return -self.zorder3d + vs[2].mean() / cz.ptp()
 
             def draw(self, renderer):
                 if self.reflen:
@@ -325,7 +325,7 @@ if mpl_enabled:
                     factor = proj_len * (args[0] +
                                          args[3]) * 0.5 * 72.0 / self.figure.dpi
 
-                    self.set_linewidths(self._linewidths_orig2 * factor)
+                    self.set_linewidths(self.linewidths_orig2 * factor)
 
                 super(Path3DCollection, self).draw(renderer)
 
diff --git a/kwant/solvers/common.py b/kwant/solvers/common.py
index a6bdecdaed1f9a693bc82ac9ced86211b9731ab6..e8655648ca1f9f661038d82a470056924c7aed50 100644
--- a/kwant/solvers/common.py
+++ b/kwant/solvers/common.py
@@ -560,11 +560,11 @@ class BlockResult(object):
         self.lead_info = lead_info
         self.out_leads = out_leads
         self.in_leads = in_leads
-        self._sizes = np.array(sizes)
-        self._in_offsets = np.zeros(len(self.in_leads) + 1, int)
-        self._in_offsets[1 :] = np.cumsum(self._sizes[self.in_leads])
-        self._out_offsets = np.zeros(len(self.out_leads) + 1, int)
-        self._out_offsets[1 :] = np.cumsum(self._sizes[self.out_leads])
+        self.sizes = np.array(sizes)
+        self.in_offsets = np.zeros(len(self.in_leads) + 1, int)
+        self.in_offsets[1 :] = np.cumsum(self.sizes[self.in_leads])
+        self.out_offsets = np.zeros(len(self.out_leads) + 1, int)
+        self.out_offsets[1 :] = np.cumsum(self.sizes[self.out_leads])
 
     def block_coords(self, lead_out, lead_in):
         """
@@ -576,16 +576,16 @@ class BlockResult(object):
         """Return a slice with the rows in the block corresponding to lead_out.
         """
         lead_out = self.out_leads.index(lead_out)
-        return slice(self._out_offsets[lead_out],
-                     self._out_offsets[lead_out + 1])
+        return slice(self.out_offsets[lead_out],
+                     self.out_offsets[lead_out + 1])
 
     def in_block_coords(self, lead_in):
         """
         Return a slice with the columns in the block corresponding to lead_in.
         """
         lead_in = self.in_leads.index(lead_in)
-        return slice(self._in_offsets[lead_in],
-                     self._in_offsets[lead_in + 1])
+        return slice(self.in_offsets[lead_in],
+                     self.in_offsets[lead_in + 1])
 
     def submatrix(self, lead_out, lead_in):
         """Return the matrix elements from lead_in to lead_out."""
diff --git a/kwant/solvers/tests/test_mumps.py b/kwant/solvers/tests/test_mumps.py
index fff9cdc1539eef56e8b44b428d621751ec12d2a5..c0f4007841e3a809521352dd1be6aeaf595ce05d 100644
--- a/kwant/solvers/tests/test_mumps.py
+++ b/kwant/solvers/tests/test_mumps.py
@@ -12,9 +12,9 @@ try:
     from kwant.solvers.mumps import \
         smatrix, greens_function, ldos, wave_function, options, reset_options
     import _test_sparse
-    _no_mumps = False
+    no_mumps = False
 except ImportError:
-    _no_mumps = True
+    no_mumps = True
 
 
 opt_list=[{},
@@ -25,7 +25,7 @@ opt_list=[{},
           {'nrhs' : 2, 'ordering' : 'amd', 'sparse_rhs' : True}]
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_output():
     for opts in opt_list:
         reset_options()
@@ -33,7 +33,7 @@ def test_output():
         _test_sparse.test_output(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_one_lead():
     for opts in opt_list:
         reset_options()
@@ -41,7 +41,7 @@ def test_one_lead():
         _test_sparse.test_one_lead(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_smatrix_shape():
     for opts in opt_list:
         reset_options()
@@ -49,14 +49,14 @@ def test_smatrix_shape():
         _test_sparse.test_smatrix_shape(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_two_equal_leads():
     for opts in opt_list:
         reset_options()
         options(**opts)
         _test_sparse.test_two_equal_leads(smatrix)
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_graph_system():
     for opts in opt_list:
         reset_options()
@@ -64,7 +64,7 @@ def test_graph_system():
         _test_sparse.test_graph_system(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_singular_graph_system():
     for opts in opt_list:
         reset_options()
@@ -72,7 +72,7 @@ def test_singular_graph_system():
         _test_sparse.test_singular_graph_system(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_tricky_singular_hopping():
     for opts in opt_list:
         reset_options()
@@ -80,7 +80,7 @@ def test_tricky_singular_hopping():
         _test_sparse.test_tricky_singular_hopping(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_selfenergy():
     for opts in opt_list:
         reset_options()
@@ -88,7 +88,7 @@ def test_selfenergy():
         _test_sparse.test_selfenergy(greens_function, smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_selfenergy_reflection():
     for opts in opt_list:
         reset_options()
@@ -96,7 +96,7 @@ def test_selfenergy_reflection():
         _test_sparse.test_selfenergy_reflection(greens_function, smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_very_singular_leads():
     for opts in opt_list:
         reset_options()
@@ -104,7 +104,7 @@ def test_very_singular_leads():
         _test_sparse.test_very_singular_leads(smatrix)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_ldos():
     for opts in opt_list:
         reset_options()
@@ -112,7 +112,7 @@ def test_ldos():
         _test_sparse.test_ldos(ldos)
 
 
-@skipif(_no_mumps)
+@skipif(no_mumps)
 def test_wavefunc_ldos_consistency():
     for opts in opt_list:
         options(**opts)