diff --git a/veerer/__init__.py b/veerer/__init__.py
index 38586db..c152ac1 100644
--- a/veerer/__init__.py
+++ b/veerer/__init__.py
@@ -31,11 +31,6 @@
 from .linear_family import VeeringTriangulationLinearFamily
 from .automaton import FlipGraph, CoreAutomaton, ReducedCoreAutomaton, GeometricAutomaton, GeometricAutomatonSubspace
 from .flip_sequence import VeeringFlipSequence
-
-from .env import sage
-if sage is not None:
-    from .flat_structure import FlatVeeringTriangulation
-    from .layout import FlatVeeringTriangulationLayout
-    from .measured_train_track import MeasuredTrainTrack
-
-del sage
+from .flat_structure import FlatVeeringTriangulation
+from .layout import FlatVeeringTriangulationLayout
+from .measured_train_track import MeasuredTrainTrack
diff --git a/veerer/automaton.py b/veerer/automaton.py
index 770d991..1dd5df7 100644
--- a/veerer/automaton.py
+++ b/veerer/automaton.py
@@ -32,9 +32,7 @@
 from .linear_family import VeeringTriangulationLinearFamily
 from .constants import RED, BLUE, PURPLE, PROPERTIES_COLOURS, colour_to_char, colour_to_string
 from .permutation import perm_invert
-
-from . import env
-from .env import sage, require_package
+from .env import CHECK
 
 
 class Automaton(object):
@@ -517,7 +515,7 @@ def set_seed(self, state):
         state = self._seed_setup(state)
 
         # TODO: check to be removed
-        if env.CHECK:
+        if CHECK:
             state._check()
 
         if self._verbosity:
@@ -604,7 +602,7 @@ def run(self, max_size=None):
             new_flip = flip_branch[-1].pop()
 
             # TODO: check to be removed
-            if env.CHECK:
+            if CHECK:
                 T._check()
 
             if self._verbosity >= 2:
@@ -670,7 +668,7 @@ def run(self, max_size=None):
             self._flip_back(flip_back_data)
 
             # TODO: check to be removed
-            if env.CHECK:
+            if CHECK:
                 T._check()
             assert T == state_branch[-1], (T, state_branch[-1])
 
@@ -726,12 +724,12 @@ def _forward_flips(self, state):
 
     def _flip(self, flip_data):
         e = flip_data
-        self._state.flip(e, check=env.CHECK)
+        self._state.flip(e, check=CHECK)
         return True, e
 
     def _flip_back(self, flip_back_data):
         e = flip_back_data
-        self._state.flip_back(e, check=env.CHECK)
+        self._state.flip_back(e, check=CHECK)
 
 
 class CoreAutomaton(Automaton):
@@ -757,13 +755,13 @@ def _forward_flips(self, state):
     def _flip(self, flip_data):
         e, col = flip_data
         old_col = self._state.colour(e)
-        self._state.flip(e, col, check=env.CHECK)
+        self._state.flip(e, col, check=CHECK)
         flip_back_data = (e, old_col)
         return self._state.edge_has_curve(e), flip_back_data
 
     def _flip_back(self, flip_back_data):
         e, old_col = flip_back_data
-        self._state.flip_back(e, old_col, check=env.CHECK)
+        self._state.flip_back(e, old_col, check=CHECK)
 
 
 class ReducedCoreAutomaton(Automaton):
@@ -784,7 +782,7 @@ def _forward_flips(self, state):
         """
         ffe = state.purple_edges()
         # TODO: check to be removed
-        if env.CHECK:
+        if CHECK:
             assert ffe == state.forward_flippable_edges()
         return [(x, col) for x in ffe for col in (BLUE, RED)]
 
@@ -793,10 +791,10 @@ def _flip(self, flip_data):
         e, col = flip_data
 
         # TODO: check to be removed
-        if env.CHECK:
+        if CHECK:
             assert T.is_forward_flippable(e)
         old_col = self._state.colour(e)
-        self._state.flip(e, col, reduced=False, check=env.CHECK)
+        self._state.flip(e, col, reduced=False, check=CHECK)
 
         if not self._state.edge_has_curve(e):
             return False, (e, old_col, ())
@@ -821,7 +819,7 @@ def _flip(self, flip_data):
                 T._colouring[T._ep[d]] = PURPLE
 
             # assertions to be removed
-            if env.CHECK:
+            if CHECK:
                 assert not T.is_forward_flippable(e)
                 assert not T.is_forward_flippable(a)
                 assert not T.is_forward_flippable(c)
@@ -837,7 +835,7 @@ def _flip(self, flip_data):
                 T._colouring[T._ep[c]] = PURPLE
 
             # assertions to be removed
-            if env.CHECK:
+            if CHECK:
                 assert not T.is_forward_flippable(e)
                 assert not T.is_forward_flippable(b)
                 assert not T.is_forward_flippable(d)
@@ -849,7 +847,7 @@ def _flip_back(self, flip_back_data):
         for ee, ccol in recolorings:
             self._state._colouring[ee] = ccol
             self._state._colouring[self._state._ep[ee]] = ccol
-        self._state.flip_back(e, old_col, check=env.CHECK)
+        self._state.flip_back(e, old_col, check=CHECK)
 
 
 class GeometricAutomaton(Automaton):
@@ -900,15 +898,15 @@ def _flip(self, flip_data):
         assert all(self._state.colour(e) == self._state.colour(edges[0]) for e in edges)
         flip_back_data = (edges, self._state.colour(edges[0]))
         for e in edges:
-            self._state.flip(e, col, check=env.CHECK)
-        if env.CHECK and not self._state.is_geometric(backend=self._backend):
+            self._state.flip(e, col, check=CHECK)
+        if CHECK and not self._state.is_geometric(backend=self._backend):
             raise RuntimeError('that was indeed possible!')
         return True, flip_back_data
 
     def _flip_back(self, flip_back_data):
         edges, old_col = flip_back_data
         for e in edges:
-            self._state.flip_back(e, old_col, check=env.CHECK)
+            self._state.flip_back(e, old_col, check=CHECK)
 
 
 class GeometricAutomatonSubspace(Automaton):
@@ -970,8 +968,8 @@ def _flip(self, flip_data):
         assert all(self._state.colour(e) == self._state.colour(edges[0]) for e in edges)
         flip_back_data = (edges, self._state.colour(edges[0]))
         for e in edges:
-            self._state.flip(e, col, check=env.CHECK)
-        if env.CHECK:
+            self._state.flip(e, col, check=CHECK)
+        if CHECK:
             self._state._check(RuntimeError)
             if not self._state.is_geometric(backend=self._backend):
                 raise RuntimeError
@@ -980,4 +978,4 @@ def _flip(self, flip_data):
     def _flip_back(self, flip_back_data):
         edges, old_col = flip_back_data
         for e in edges:
-            self._state.flip_back(e, old_col, check=env.CHECK)
+            self._state.flip_back(e, old_col, check=CHECK)
diff --git a/veerer/env.py b/veerer/env.py
index 755cb9c..3c2e305 100644
--- a/veerer/env.py
+++ b/veerer/env.py
@@ -1,10 +1,5 @@
 r"""
-Various environment modules.
-
-TESTS::
-
-    sage: from veerer.env import sage, flipper, surface_dynamics, ppl  # random output due to deprecation warnings from realalg
-
+Tests for optional packages used by veerer.
 """
 # ****************************************************************************
 #  This file is part of veerer
@@ -28,61 +23,20 @@
 
 CHECK = False
 
-import sage.all
-import sage
-import sage.misc.prandom as random
-from sage.structure.richcmp import op_LT, op_LE, op_EQ, op_NE, op_GT, op_GE, rich_to_bool
-
-
-try:
-    import surface_dynamics
-except ImportError:
-    surface_dynamics = None
-
-
-try:
-    import flipper
-except ImportError:
-    flipper = None
-
-try:
-    import curver
-except ImportError:
-    curver = None
-
-try:
-    import ppl
-except ImportError:
-    ppl = None
-
-
-try:
-    import PyNormaliz
-except ImportError:
-    PyNormaliz = None
-
-error_msg = {
-    'curver': 'the function {} can only be called when the package curver is installed.',
-
-    'surface_dynamics': 'the function {} only works when the package surface_dynamics is installed. See https://pypi.org/project/surface_dynamics/ for instructions.',
-
-    'flipper': 'the function {} only works when the package flipper is installed. See https://pypi.org/project/flipper/ for instructions',
-
-    'ppl': 'the function {} only works when the package pplpy is installed. See https://pypi.org/project/pplpy/ for instructions.',
-
-    'PyNormaliz': 'the function {} only works when the package PyNormaliz is installed.'
-    }
-
-missing_mods = {
-    'curver': curver is None,
-    'flipper': flipper is None,
-    'ppl': ppl is None,
-    'surface_dynamics': surface_dynamics is None,
-    'PyNormaliz': PyNormaliz is None
-    }
-
-# TODO: use the traceback to find out who called this function!
-# https://docs.python.org/2/library/traceback.html#traceback-examples
-def require_package(mod_name, caller):
-    if missing_mods[mod_name]:
-        raise ValueError(error_msg[mod_name].format(caller))
+from sage.features import PythonModule
+
+pyflatsurf_feature = PythonModule(
+    "pyflatsurf", url="https://github.com/flatsurf/flatsurf/#install-with-conda"
+)
+surface_dynamics_feature = PythonModule(
+    "surface-dynamics", url="https://flatsurf.github.io/surface-dynamics/#installation"
+) 
+flipper_feature = PythonModule(
+    "flipper", url="https://flipper.readthedocs.io/en/latest/user/install.html"
+)
+curver_feature = PythonModule(
+    "curver", url="https://curver.readthedocs.io/en/master/user/install.html"
+)
+pynormaliz_feature = PythonModule(
+    "PyNormaliz", url="https://github.com/Normaliz/PyNormaliz"
+)
diff --git a/veerer/features.py b/veerer/features.py
new file mode 100644
index 0000000..9e09410
--- /dev/null
+++ b/veerer/features.py
@@ -0,0 +1,21 @@
+r"""
+Tests for optional packages used by veerer.
+"""
+
+from sage.features import PythonModule
+
+flatsurf_feature = PythonModule(
+    "flatsurf", url="https://flatsurf.github.io/sage-flatsurf/#installation"
+)
+surface_dynamics_feature = PythonModule(
+    "surface_dynamics", url="https://flatsurf.github.io/surface-dynamics/#installation"
+)
+flipper_feature = PythonModule(
+    "veerer", url="https://flatsurf.github.io/veerer/installation.html"
+)
+curver_feature = PythonModule(
+    "curver", url="https://curver.readthedocs.io/en/master/user/install.html"
+)
+pynormaliz_feature = PythonModule(
+    "PyNormaliz", url="https://github.com/Normaliz/PyNormaliz"
+)
diff --git a/veerer/flip_sequence.py b/veerer/flip_sequence.py
index a2d75a2..bdb0d2e 100644
--- a/veerer/flip_sequence.py
+++ b/veerer/flip_sequence.py
@@ -59,7 +59,7 @@
 from .constants import colour_from_char, colour_to_char, RED, BLUE, PURPLE, GREEN, HORIZONTAL, VERTICAL
 from .permutation import perm_init, perm_check, perm_id, perm_is_one, perm_preimage, perm_invert, perm_cycle_string, perm_compose, perm_pow, perm_conjugate
 from .veering_triangulation import VeeringTriangulation
-from .env import require_package, sage, ppl
+
 
 def flip_sequence_to_string(sequence):
     return " ".join("%d%s" % (e, colour_to_char(col)) for e,col in sequence)
diff --git a/veerer/linear_family.py b/veerer/linear_family.py
index 41910b0..f3ad178 100644
--- a/veerer/linear_family.py
+++ b/veerer/linear_family.py
@@ -28,33 +28,21 @@
 import numbers
 from random import choice, shuffle
 
-from .env import sage, ppl
+from sage.structure.element import get_coercion_model, Matrix
+from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+from sage.matrix.constructor import matrix
+from sage.geometry.polyhedron.constructor import Polyhedron
+from sage.arith.misc import gcd
+from sage.categories.number_fields import NumberFields
+
 from .constants import VERTICAL, HORIZONTAL, BLUE, RED
 from .permutation import perm_cycle_string, perm_cycles, perm_check, perm_conjugate, perm_on_list
 from .polyhedron import LinearExpressions, ConstraintSystem
 from .veering_triangulation import VeeringTriangulation
 
-
-if sage is not None:
-    from sage.structure.element import get_coercion_model
-    from sage.rings.integer_ring import ZZ
-    from sage.rings.rational_field import QQ
-    from sage.matrix.constructor import matrix
-    from sage.geometry.polyhedron.constructor import Polyhedron
-    from sage.arith.misc import gcd
-    from sage.categories.number_fields import NumberFields
-
-    cm = get_coercion_model()
-    _NumberFields = NumberFields()
-else:
-    matrix = None
-    Polyhedron = None
-    gcd = None
-    ZZ = None
-    QQ = None
-
-    cm = None
-    _NumberFields = None
+cm = get_coercion_model()
+_NumberFields = NumberFields()
 
 
 def subspace_are_equal(subspace1, subspace2, check=True):
@@ -210,7 +198,7 @@ def __init__(self, *args, mutable=False, check=True):
             t, colouring, subspace = args
         VeeringTriangulation.__init__(self, t, colouring, mutable=mutable, check=False)
 
-        if not isinstance(subspace, sage.structure.element.Matrix):
+        if not isinstance(subspace, Matrix):
             subspace = matrix(subspace)
 
         self._subspace = subspace
diff --git a/veerer/misc.py b/veerer/misc.py
index df14717..a44d0d2 100644
--- a/veerer/misc.py
+++ b/veerer/misc.py
@@ -21,11 +21,13 @@
 #  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 # ****************************************************************************
 
-from .env import ppl
+import ppl
+
 
 def det2(u, v):
     return u[0]*v[1] - u[1]*v[0]
 
+
 def flipper_edge(T, e):
     r"""
     EXAMPLES::
@@ -40,10 +42,12 @@ def flipper_edge(T, e):
     e = e.label
     return n * (e < 0) + e
 
+
 def flipper_edge_perm(n):
     from array import array
     return array('i', range(n-1,-1,-1))
 
+
 def flipper_face_edge_perms(T):
     r"""
     Return a pair ``(face permutation, edge permutation)`` from a flipper triangulation.
@@ -55,6 +59,7 @@ def flipper_face_edge_perms(T):
     triangles = [(flipper_edge(T, e), flipper_edge(T, f), flipper_edge(T, g)) for e,f,g in T]
     return perm_init(triangles), flipper_edge_perm(n)
 
+
 def flipper_isometry_to_perm(isom, ep, inv=False):
     r"""
     Question: how do we create an isometry in flipper?
@@ -75,6 +80,7 @@ def flipper_isometry_to_perm(isom, ep, inv=False):
         p[i] = j
     return p
 
+
 def flipper_nf_to_sage(K, name='a'):
     r"""
     Convert a flipper number field to Sage.
@@ -97,6 +103,7 @@ def flipper_nf_to_sage(K, name='a'):
     s = AA.polynomial_root(p, RIF(l,u))
     return NumberField(p, name, embedding=s)
 
+
 def flipper_nf_element_to_sage(x, K=None):
     r"""
     Convert a flipper nf element to Sage.
@@ -120,6 +127,7 @@ def flipper_nf_element_to_sage(x, K=None):
         coeffs.extend([0] * (K.degree() - len(coeffs)))
     return K(list(map(QQ, coeffs)))
 
+
 def rays_to_ppl_cone(rays):
     gs = ppl.Generator_System()
     gs.insert(ppl.point())
diff --git a/veerer/polyhedron/cone.py b/veerer/polyhedron/cone.py
index d72b2f0..09e1090 100644
--- a/veerer/polyhedron/cone.py
+++ b/veerer/polyhedron/cone.py
@@ -86,38 +86,32 @@
 #  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 # ****************************************************************************
 
-from ..env import sage, ppl, PyNormaliz
-if sage is not None:
-    from sage.categories.number_fields import NumberFields
+from fractions import Fraction
+
+from sage.categories.number_fields import NumberFields
+from sage.rings.integer_ring import ZZ
+from sage.rings.rational_field import QQ
+from sage.rings.real_arb import RealBallField
+from sage.rings.real_double import RDF
+from sage.rings.qqbar import AA, number_field_elements_from_algebraics
 
-    from sage.rings.integer_ring import ZZ
-    from sage.rings.rational_field import QQ
-    from sage.rings.real_arb import RealBallField
-    from sage.rings.real_double import RDF
-    from sage.rings.qqbar import AA, number_field_elements_from_algebraics
+from ..features import pynormaliz_feature
 
-    RBF = RealBallField(64)
-    _NumberFields = NumberFields()
 
+RBF = RealBallField(64)
+_NumberFields = NumberFields()
 
 _backends = {}
-from fractions import Fraction
-if ppl is not None:
-    _backends[int] = _backends[Fraction] = 'ppl'
-elif PyNormaliz is not None:
-    _backends[int] = _backends[Fraction] = 'normaliz-QQ'
-if sage is not None:
-    from sage.rings.integer_ring import ZZ
-    from sage.rings.rational_field import QQ
-    _backends[ZZ] = _backends[QQ] = 'ppl'
+
+_backends[int] = _backends[Fraction] = 'ppl'
+_backends[ZZ] = _backends[QQ] = 'ppl'
 def default_backend(base_ring):
     try:
         return _backends[base_ring]
     except KeyError:
         pass
-    if sage is not None:
-        if base_ring in _NumberFields:
-            return 'normaliz-NF' if PyNormaliz is not None else 'sage'
+    if base_ring in _NumberFields:
+        return 'normaliz-NF' if pynormaliz_feature.is_present() else 'sage'
     raise ValueError('no suitable polytope backend for base ring {}'.format(base_ring))
 
 
diff --git a/veerer/polyhedron/linear_expression.py b/veerer/polyhedron/linear_expression.py
index 348b5a7..75f8b26 100644
--- a/veerer/polyhedron/linear_expression.py
+++ b/veerer/polyhedron/linear_expression.py
@@ -24,18 +24,15 @@
 #  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
 # ****************************************************************************
 
-from ..env import require_package, sage, ppl
+import ppl
 
 from sage.categories.modules import Modules
-
 from sage.structure.unique_representation import UniqueRepresentation
 from sage.structure.element import ModuleElement, Vector, parent, get_coercion_model
 from sage.structure.parent import Parent
 from sage.structure.richcmp import op_LE, op_LT, op_EQ, op_NE, op_GT, op_GE
-
 from sage.rings.integer_ring import ZZ
 from sage.rings.rational_field import QQ
-
 from sage.arith.functions import lcm
 from sage.arith.functions import LCM_list
 
@@ -227,7 +224,6 @@ def ppl(self):
             sage: (3 * L.variable(2) - 7/2 * L.variable(5) + 1/3).ppl()
             18*x2-21*x5+2
         """
-        require_package('ppl', 'ppl')
         from gmpy2 import mpz
         lin = self.integral()
         # TODO: the line below is too costly : it accounts for 80% of the
@@ -557,7 +553,6 @@ def cone(self, backend='sage'):
             backend = default_backend(self.base_ring())
 
         if backend == 'ppl':
-            require_package('ppl', 'cone')
             from .cone import Cone_ppl
             return Cone_ppl(QQ, ppl.C_Polyhedron(self.ppl()))
         elif backend == 'sage':
@@ -565,12 +560,16 @@ def cone(self, backend='sage'):
             ieqs, eqns = self.ieqs_eqns()
             return Cone_sage._new(ieqs, eqns)
         elif backend == 'normaliz-QQ':
-            require_package('PyNormaliz', 'cone')
+            from .features import pynormaliz_feature
+            pynormaliz_feature.require()
+
             from .cone import Cone_normalizQQ
             ieqs, eqns = self.ieqs_eqns(homogeneous=True)
             return Cone_normalizQQ._new(ieqs, eqns)
         elif backend == 'normaliz-NF':
-            require_package('PyNormaliz', 'cone')
+            from .features import pynormaliz_feature
+            pynormaliz_feature.require()
+
             from .cone import Cone_normalizNF
             ieqs, eqns = self.ieqs_eqns(homogeneous=True)
             return Cone_normalizNF._new(ieqs, eqns)
diff --git a/veerer/triangulation.py b/veerer/triangulation.py
index 001c4ae..a278dab 100644
--- a/veerer/triangulation.py
+++ b/veerer/triangulation.py
@@ -26,12 +26,13 @@
 import numbers
 from array import array
 
+from sage.structure.richcmp import op_LT, op_LE, op_EQ, op_NE, op_GT, op_GE, rich_to_bool
+
 from .permutation import (perm_init, perm_check, perm_cycles, perm_dense_cycles,
                           perm_invert, perm_conjugate, perm_cycle_string, perm_cycles_lengths,
                           perm_num_cycles, str_to_cycles, perm_compose, perm_from_base64_str,
                           uint_base64_str, uint_from_base64_str, perm_base64_str,
                           perms_are_transitive, triangulation_relabelling_from)
-from .env import require_package, flipper, curver, sage, rich_to_bool, op_LE, op_LT, op_EQ, op_NE, op_GT, op_GE
 
 
 def face_edge_perms_init(data):
@@ -591,7 +592,10 @@ def to_flipper(self):
             sage: V.to_flipper()  # optional - flipper
             [(~2, ~0, ~1), (0, 1, 2)]
         """
-        require_package('flipper', 'to_flipper')
+        from .features import flipper_feature
+        flipper_feature.require()
+
+        import flipper
         ep = self._ep
         F = []
         for f in self.faces():
@@ -625,7 +629,10 @@ def to_curver(self):
             sage: V.to_curver()  # optional - curver
             [(~2, ~0, ~1), (0, 1, 2)]
         """
-        require_package('curver', 'to_curver')
+        from .features import curver_feature
+        curver_feature.require()
+
+        import curver
         ep = self._ep
         F = []
         for f in self.faces():
diff --git a/veerer/veering_quadrangulation.py b/veerer/veering_quadrangulation.py
index 0d80afb..9aa2908 100644
--- a/veerer/veering_quadrangulation.py
+++ b/veerer/veering_quadrangulation.py
@@ -48,7 +48,6 @@
 from sage.matrix.special import identity_matrix, zero_matrix, block_matrix
 from sage.libs.gap.libgap import libgap
 
-from veerer.env import require_package
 from veerer.permutation import *
 
 RIGHT = 1
@@ -157,7 +156,9 @@ def __repr__(self):
                  perm_cycle_string(self._pl, n=self._n))
 
     def to_origami(self):
-        require_package('surface_dynamics', 'to_origami')
+        from .features import surface_dynamics_feature
+        surface_dynamics_feature.require()
+
         from surface_dynamics import Origami
         return Origami(list(self._pr), list(self._pl), as_tuple=True, check=False)
 
diff --git a/veerer/veering_triangulation.py b/veerer/veering_triangulation.py
index febb754..893f044 100644
--- a/veerer/veering_triangulation.py
+++ b/veerer/veering_triangulation.py
@@ -32,6 +32,9 @@
 import numbers
 from random import choice, shuffle
 from array import array
+import ppl
+
+from sage.structure.richcmp import op_LT, op_LE, op_EQ, op_NE, op_GT, op_GE, rich_to_bool
 
 from .constants import *
 from .permutation import *
@@ -40,8 +43,6 @@
 from .polyhedron import LinearExpressions, ConstraintSystem
 from .polyhedron.linear_algebra import linear_form_project, linear_form_normalize
 
-from .env import curver, sage, surface_dynamics, ppl, flipper, random, require_package, rich_to_bool, op_LE, op_LT, op_EQ, op_NE, op_GT, op_GE, CHECK
-
 
 class VeeringTriangulation(Triangulation):
     r"""
@@ -434,7 +435,9 @@ def from_square_tiled(cls, s, col=RED, mutable=False, check=True):
             sage: T.stratum()                                     # optional - surface_dynamics
             H_3(4)
         """
-        require_package('surface_dynamics', 'from_square_tiled')
+        from .features import surface_dynamics_feature
+        surface_dynamics_feature.require()
+
         from surface_dynamics.flat_surfaces.origamis.origami_dense import Origami_dense_pyx
 
         if col not in [BLUE, RED]:
@@ -516,7 +519,8 @@ def from_stratum(cls, c, folded_edges=False, mutable=False, check=True):
             sage: CT.stratum()                                                # optional - surface_dynamics
             H_4(6)
         """
-        require_package('surface_dynamics', 'from_stratum')
+        from .features import surface_dynamics_feature
+        surface_dynamics_feature.require()
 
         # TODO: for now there is no account of possible folded edges
         # in the cylinder diagram. This has to be changed in
@@ -1124,7 +1128,8 @@ def stratum(self):
             sage: t.stratum()                                                    # optional - surface_dynamics
             Q_2(1^4)
         """
-        require_package('surface_dynamics', 'stratum')
+        from .features import surface_dynamics_feature
+        surface_dynamics_feature.require()
 
         A = self.angles()
         if any(a%2 for a in A) or not self.is_abelian():
@@ -2945,8 +2950,6 @@ def train_track_min_solution(self, slope=VERTICAL, allow_degenerations=False):
             sage: T.train_track_min_solution(HORIZONTAL, allow_degenerations=True)
             point(1/1, 0/1, 1/1)
         """
-        require_package('ppl', 'train_track_min_solution')
-
         n = self.num_edges()
         M = ppl.MIP_Problem(n)