diff --git a/qsymm/__init__.py b/qsymm/__init__.py
index 43bd01ed7969696796752962450c4c3d3e4a160e..65b78de8850f651197e92adfbab6e43dce78589a 100644
--- a/qsymm/__init__.py
+++ b/qsymm/__init__.py
@@ -13,5 +13,3 @@ from .hamiltonian_generator import continuum_hamiltonian, continuum_pairing, dis
 from .symmetry_finder import symmetries, discrete_symmetries, conserved_quantities, \
                              continuous_symmetries, bravais_point_group
 from .kwant_continuum import sympify
-
-del _version
diff --git a/qsymm/groups.py b/qsymm/groups.py
index 5160d904a9ed31d3f20294d24b68cfc496c0d3f8..34eafd14d59962b25c6ee6f19335cbefc2ac51ab 100644
--- a/qsymm/groups.py
+++ b/qsymm/groups.py
@@ -13,7 +13,6 @@ import sympy
 from sympy.matrices.matrices import MatrixBase
 
 from .linalg import prop_to_id, _inv_int, allclose
-from .model import Model
 from .kwant_continuum import sympify
 
 
@@ -755,7 +754,7 @@ def cubic(tr=True, ph=True, generators=False, spin=None):
         raise ValueError('If `ph` is True, `spin` may not be provided, as it is not '
                          'possible to deduce the unitary representation of particle-hole symmetry '
                          'from spin alone. In this case construct the particle-hole operator manually.')
-    I = inversion(3, U=(None if spin is None else spin_rotation(np.zeros(3), spin)))
+    I = inversion(3, U=(None if spin is None else spin_rotation(np.zeros(3), spin)))  # Noqa: E741
     C4 = rotation(1/4, [1, 0, 0], spin=spin)
     C3 = rotation(1/3, [1, 1, 1], spin=spin)
     cubic_gens = {I, C4, C3}
@@ -829,7 +828,7 @@ def hexagonal(dim=2, tr=True, ph=True, generators=False, sympy_R=True, spin=None
             C6 = rotation(1/6, spin=spin)
         gens = {Mx, C6}
     elif dim == 3:
-        I = inversion(3, U=(None if spin is None else spin_rotation(np.zeros(3), spin)))
+        I = inversion(3, U=(None if spin is None else spin_rotation(np.zeros(3), spin)))  # Noqa: E741
         C2x = rotation(1/2, [1, 0, 0], spin=spin)
         if sympy_R:
             C6 = PointGroupElement(sympy.ImmutableMatrix(
@@ -1037,7 +1036,7 @@ def pretty_print_cgg(g, latex=False):
         L = None
 
     if R is not None and R.shape[0] == 3:
-        n = np.array([np.trace(l @ R) for l in L_matrices()]).real
+        n = np.array([np.trace(l @ R) for l in L_matrices()]).real  # Noqa: E741
         n /= la.norm(n)
         n = _round_axis(n)
         if latex:
@@ -1195,7 +1194,7 @@ def spin_rotation(n, s, roundint=False):
     return U
 
 
-def L_matrices(d=3, l=1):
+def L_matrices(d=3, l=1):  # Noqa: E741
     """Construct real space rotation generator matrices in d=2 or 3 dimensions.
     Can also be used to get angular momentum operators for real atomic orbitals
     in 3 dimensions, for p-orbitals use `l=1`, for d-orbitals `l=2`. The basis
diff --git a/qsymm/hamiltonian_generator.py b/qsymm/hamiltonian_generator.py
index 031547698c7f353ede272f80e6db34ba3597e881..ebda7e9e28007254d3cf164dbb8741e08aafe2c2 100644
--- a/qsymm/hamiltonian_generator.py
+++ b/qsymm/hamiltonian_generator.py
@@ -6,7 +6,7 @@ import scipy.linalg as la
 from copy import deepcopy
 import tinyarray as ta
 
-from .linalg import matrix_basis, nullspace, sparse_basis, family_to_vectors, rref, allclose
+from .linalg import matrix_basis, nullspace, family_to_vectors, rref, allclose
 from .model import Model, BlochModel, BlochCoeff, _commutative_momenta, e, I
 from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group
 from . import kwant_continuum
@@ -55,7 +55,7 @@ def continuum_hamiltonian(symmetries, dim, total_power, all_powers=None,
         all the symmetries by construction.
     """
 
-    if type(total_power) is int:
+    if isinstance(total_power, int):
         max_power = total_power
         total_power = range(max_power + 1)
 
@@ -132,11 +132,11 @@ def continuum_pairing(symmetries, dim, total_power, all_powers=None, momenta=_co
         satisfies the symmetries specified. Each Model object satisfies
         all the symmetries by construction.
     """
-    if type(total_power) is int:
+    if isinstance(total_power, int):
         max_power = total_power
         total_power = range(max_power + 1)
 
-    if not ph_square in (-1, 1):
+    if ph_square not in (-1, 1):
         raise ValueError('Particle-hole operator must square to +1 or -1.')
     if phases is None:
         phases = np.ones(len(symmetries))
@@ -218,13 +218,13 @@ def continuum_variables(dim, total_power, all_powers=None, momenta=_commutative_
     all_powers = all_powers[:dim]
 
     for i, power in enumerate(all_powers):
-        if type(power) is int:
+        if isinstance(power, int):
             all_powers[i] = range(power + 1)
 
     if dim == 0:
         return [kwant_continuum.sympify(1)]
 
-    if all([type(i) is int for i in momenta]):
+    if all(isinstance(i, int) for i in momenta):
         momenta = [_commutative_momenta[i] for i in momenta]
     else:
         momenta = [kwant_continuum.make_commutative(k, k)
diff --git a/qsymm/kwant_linalg_lll.py b/qsymm/kwant_linalg_lll.py
index 0d81ddf2d14451486d2c7970ce0aef83420bb3cc..2f10c1e345a2ad3ce255374cfd493d427b92efb8 100644
--- a/qsymm/kwant_linalg_lll.py
+++ b/qsymm/kwant_linalg_lll.py
@@ -206,7 +206,7 @@ def cvp(vec, basis, n=1, group_by_length=False, rtol=1e-09):
     # size as the diameter of the largest sphere inside.
     rad = 0.5 / np.sqrt(np.max(np.diag(la.inv(bbt))))
 
-    l = 1
+    l = 1  # Noqa: E741
     while True:
         # Make all the lattice points in and on the edges of a parallelepiped
         # of `2*l` size around `center_coords`.
@@ -214,7 +214,7 @@ def cvp(vec, basis, n=1, group_by_length=False, rtol=1e-09):
         points = points.reshape(basis.shape[0], -1).T
         # If there are less than `n` points, make more.
         if len(points) < n:
-            l += 1
+            l += 1  # Noqa: E741
             continue
         point_coords = points @ basis
         point_coords = point_coords - vec.T
@@ -227,7 +227,7 @@ def cvp(vec, basis, n=1, group_by_length=False, rtol=1e-09):
             group_boundaries = np.where(np.diff(distances) > rtol * rad)[0]
             if len(group_boundaries) + 1 < n:
                 # Make more points if there are less than `n` groups.
-                l += 1
+                l += 1  # Noqa: E741
                 continue
             elif len(group_boundaries) == n - 1:
                 # If there are exactly `n` groups, choose largest distance.
@@ -255,7 +255,7 @@ def cvp(vec, basis, n=1, group_by_length=False, rtol=1e-09):
         else:
             # Otherwise there may be smaller vectors we haven't found,
             # increase `l`.
-            l += 1
+            l += 1  # Noqa: E741
             continue
 
 
diff --git a/qsymm/model.py b/qsymm/model.py
index f0d7522c8ed6551870fc92bd25f2e23c425c0441..fdb6b1d4f144faa8eea649aa5280d964e8e15c4b 100644
--- a/qsymm/model.py
+++ b/qsymm/model.py
@@ -717,7 +717,8 @@ class Model(UserDict):
         if self.data == {}:
             return result
         # Write it explicitely so it works with sparse arrays
-        norm = lambda mat: np.sqrt(np.sum(np.abs(mat)**2))
+        def norm(mat):
+            return np.sqrt(np.sum(np.abs(mat)**2))
         max_norm = np.max([norm(val) for val in self.values()])
         tol = max(atol, max_norm * rtol)
         result.data = {key: copy(val) for key, val in self.items() if not norm(val) < tol}
diff --git a/qsymm/symmetry_finder.py b/qsymm/symmetry_finder.py
index 8e3a0e2b8e6700e43799c7d952ccc7c60bfee7d5..3dfc11e8782e40f8a05d8b11a2d92c80ee887628 100644
--- a/qsymm/symmetry_finder.py
+++ b/qsymm/symmetry_finder.py
@@ -9,13 +9,13 @@ import scipy.sparse
 from .linalg import matrix_basis, nullspace, split_list, simult_diag, commutator, \
                     prop_to_id, sparse_basis, mtm, family_to_vectors, solve_mat_eqn, \
                     allclose
-from .model import Model, BlochModel, BlochCoeff
+from .model import BlochModel, BlochCoeff
 from .groups import PointGroupElement, ContinuousGroupGenerator, generate_group, \
-                    set_multiply, L_matrices, spin_rotation, time_reversal, \
+                    set_multiply, time_reversal, \
                     particle_hole, chiral, inversion, rotation, mirror, PrettyList
 
 from . import kwant_continuum
-from .kwant_linalg_lll import lll, cvp, voronoi
+from .kwant_linalg_lll import lll, voronoi
 
 ### Top level function for symmetry finding
 
@@ -351,7 +351,7 @@ def conserved_quantities(Ps, prettify=False, num_digits=3):
         # First block does not have identity, so the identity is not
         # among the conserved quantities
         bas = matrix_basis(P.shape[0], traceless=(i==0))
-        for l in bas:
+        for l in bas:  # Noqa: E741
             # construct conserved L that acts with l on all the subblocks
             # of the original space at the same time (sum over j)
             # 'aij,ab,bkj->ik'
@@ -692,7 +692,8 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li
     if dim <= 1:
         # There is no continuous rotation in 0 and 1 D
         return []
-    Rs = lambda: matrix_basis(dim, antihermitian=True, real=True)
+    def Rs():
+        return matrix_basis(dim, antihermitian=True, real=True)
     # length of Rs
     NR = dim * (dim - 1) / 2
     blockdims = [list(rham.values())[0].shape[0] for rham in reduced_hamiltonians]
@@ -745,7 +746,7 @@ def continuous_symmetries(model, Ps=None, prettify=True, num_digits=8, sparse_li
             if blockdims[i] > 1:
                 Ls = matrix_basis(blockdims[i], traceless=True)
                 blockind = int(NR + np.sum([d**2-1 for d in blockdims[:i]]))
-                l = sum(l * v[blockind + j] for j, l in enumerate(Ls))
+                l = sum(l * v[blockind + j] for j, l in enumerate(Ls))  # Noqa: E741
                 L += np.einsum('aij,jl,akl->ik', Ps[i], l, Ps[i].conj())
         g = ContinuousGroupGenerator(R, L)
         symmetries.append(g)
diff --git a/qsymm/tests/test_groups.py b/qsymm/tests/test_groups.py
index 55032f986d9ac8d8a6f29d169820833f6ca53553..427ef4dbaef9c80f48a46be0bd32997247e3476f 100644
--- a/qsymm/tests/test_groups.py
+++ b/qsymm/tests/test_groups.py
@@ -1,9 +1,8 @@
 import pytest
-import sympy
 import numpy as np
 import tinyarray as ta
 
-from ..groups import PointGroupElement, Model, time_reversal, chiral, rotation
+from ..groups import PointGroupElement
 
 
 @pytest.mark.parametrize(
diff --git a/qsymm/tests/test_hamiltonian_generator.py b/qsymm/tests/test_hamiltonian_generator.py
index ac1503fb615b94621510fab3ae766336288c7e1a..3563fd58322a1d714e1158f3fc91f7ddd0427c86 100644
--- a/qsymm/tests/test_hamiltonian_generator.py
+++ b/qsymm/tests/test_hamiltonian_generator.py
@@ -1,4 +1,3 @@
-import pytest
 import sympy
 import numpy as np
 import scipy.linalg as la
@@ -7,7 +6,7 @@ from .. import kwant_rmt
 from ..hamiltonian_generator import continuum_hamiltonian, check_symmetry, \
      bloch_family, make_basis_pretty, constrain_family, continuum_variables, \
      continuum_pairing, remove_duplicates, subtract_family, display_family
-from ..groups import PointGroupElement, Model, time_reversal, chiral, rotation
+from ..groups import PointGroupElement, time_reversal, chiral, rotation
 from ..model import _commutative_momenta, Model, BlochModel
 from ..linalg import nullspace, family_to_vectors
 
diff --git a/qsymm/tests/test_model.py b/qsymm/tests/test_model.py
index a58e4130d12c1fcace7e7da6d3fd8121d33760b3..61875e57619605f0d58b146b4ed570431d78e7f1 100644
--- a/qsymm/tests/test_model.py
+++ b/qsymm/tests/test_model.py
@@ -1,5 +1,3 @@
-import pytest
-import warnings
 import sympy
 from sympy.core.numbers import One
 from scipy.sparse import csr_matrix
diff --git a/qsymm/tests/test_mutual.py b/qsymm/tests/test_mutual.py
index 5e42ea23ded4e2a6542b321c1da7f9b6d7f6b169..76402255e4ef4689f4399c562636cd07239c6c50 100644
--- a/qsymm/tests/test_mutual.py
+++ b/qsymm/tests/test_mutual.py
@@ -1,4 +1,3 @@
-import pytest
 import numpy as np
 import scipy.linalg as la
 from copy import deepcopy
@@ -71,7 +70,7 @@ def test_mutual_continuum():
     # More realistic symmetry action
     TR = PointGroupElement(np.eye(dim), True, False, np.kron(sigma[1], np.eye(2)))
     PH = PointGroupElement(np.eye(dim), True, True, np.kron(np.eye(2), sigma[2]))
-    I = PointGroupElement(-np.eye(dim), False, False, np.eye(4))
+    I = PointGroupElement(-np.eye(dim), False, False, np.eye(4))  # Noqa: E741
     gens = {TR, PH, I}
     group = generate_group(gens)
     groupnoU = deepcopy(group)
diff --git a/qsymm/tests/test_symmetry_finder.py b/qsymm/tests/test_symmetry_finder.py
index 36aa2a01ab0d1dd84ce85ed277f8037da0bb58d2..9ee8127b5c9ec0ceae854c191e71e9f266c607ee 100644
--- a/qsymm/tests/test_symmetry_finder.py
+++ b/qsymm/tests/test_symmetry_finder.py
@@ -1,11 +1,15 @@
 import itertools as it
-import pytest
 import sympy
+import numpy as np
+from scipy import linalg as la
 
 from .. import kwant_rmt
-from ..symmetry_finder import *
-from ..symmetry_finder import _reduced_model, _reduce_hamiltonian, bravais_point_group
-from ..linalg import *
+from ..symmetry_finder import discrete_symmetries, symmetry_adapted_sun, solve_mat_eqn,\
+    time_reversal, particle_hole, chiral, generate_group, continuous_symmetries,\
+    _reduce_hamiltonian, bravais_point_group
+from ..model import Model
+from ..groups import PointGroupElement
+from ..linalg import mtm, simult_diag
 from .. import kwant_continuum
 
 sigma = np.array([[[1, 0], [0, 1]], [[0, 1], [ 1, 0]], [[0, -1j], [1j, 0]], [[1, 0], [0, -1]]])
@@ -455,8 +459,7 @@ def test_continuum():
 
     # Cubic point group
     eye = np.array(np.eye(3), int)
-    E = PointGroupElement(eye, False, False, None)
-    I = PointGroupElement(-eye, False, False, None)
+    I = PointGroupElement(-eye, False, False, None)  # Noqa: E741
     C4 = PointGroupElement(np.array([[1, 0, 0],
                                     [0, 0, 1],
                                     [0, -1, 0]], int), False, False, None)
@@ -524,7 +527,6 @@ def test_bloch():
     # Simple tests for Bloch models
 
     # Hexagonal point group
-    eyesym = sympy.ImmutableMatrix(sympy.eye(2))
     Mx = PointGroupElement(sympy.ImmutableMatrix([[-1, 0],
                                                   [0, 1]]),
                             False, False, None)
diff --git a/qsymm/tests/test_util.py b/qsymm/tests/test_util.py
index 87b1ef1f02d4a6d36c9d3952d4e68ddccfc23792..b8340a9f625e956cd19674fc74bf8a8e4716d61e 100644
--- a/qsymm/tests/test_util.py
+++ b/qsymm/tests/test_util.py
@@ -1,12 +1,11 @@
 import pytest
-import warnings
 import sympy
 import itertools
 import numpy as np
 import scipy.linalg as la
 
 from pytest import raises
-from ..linalg import matrix_basis, sparse_basis
+from ..linalg import sparse_basis
 from ..groups import PointGroupElement
 
 
@@ -49,4 +48,4 @@ def test_spatial_types():
     assert C6s == C6f
     # Mixing sympy with other types raises an error
     with raises(ValueError):
-        S = C6s * C6f
+        C6s * C6f